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Chapter  1 

SUMMARY 


Advances  in  Air  Force  systems  require  commensurate  advancement  of  the  understanding  of  the  pacing  physics 
and  technology  challenges.  Some  of  these  challenges  involve  the  propagation  of  wave  phenomena  through 
atmospheric  media.  Electro-magnetic  responses  are  key  to  the  performance  of  sensors  and  stealth  technology. 
Acoustic  responses  will  drive  observability  and  component  fatigue.  The  Air  Force  needs  an  advanced  set  of 
computational  tools  for  computational  discovery  and  numerical  analysis  of  wave  propagation  as  it  relates  to 
flight  vehicle  performance. 

For  the  past  14  years,  scientists  in  the  Computational  Sciences  Branch  of  the  Air  Force  Research  Lab¬ 
oratory  have  systematically  pursued  a  research  program  to  develop,  validate  and  demonstrate  numerical 
approaches  for  computational  electromagnetics/acoustics  in  the  time  domain.  These  efforts  span  research  in 
high  resolution  techniques  and  methods  for  applying  them  to  Air  Force  relevant  challenges,  as  well  as  ap¬ 
plication  to  both  canonical  and  representative  configurations  and  regimes.  This  technical  report  documents 
the  development  and  demonstration  of  high-order  time-domain  computational  approaches  for  electromag¬ 
netics,  acoustics  and  aero-optic  applications.  The  methods  and  procedures  for  the  numerics  common  to  all 
applications  are  detailed  in  Chapter  3. 

Chapter  4  of  the  report  covers  details  of  the  high-orcler  methodologies  for  accurate  simulation  of  wave 
propagation  phenomena  in  general  configurations  and  though  dispersive/non-homogeneous  media.  In  order 
to  accurately  represent  the  wave  propagation,  a  6th-order  field  solver,  FDL3DI ,  was  developed  that  could 
retain  the  essential  features  of  the  wavefront.  The  impact  of  mesh  stretching  in  the  generation  of  high- 
frequency  spurious  modes  is  examined  and  the  need  for  a  discriminating  higher-order  filter  procedure  is 
established  and  resolved.  The  incorporation  of  these  filtering  techniques  also  permits  a  robust  treatment  of 
outflow  radiation  condition  by  taking  advantage  of  energy  transfer  to  high-frequencies  caused  by  rapid  mesh 
stretching.  For  conditions  on  the  scatterer,  higher-order  one-sided  filter  treatments  are  shown  to  be  superior 
in  terms  of  accuracy  and  stability  compared  to  standard  explicit  variations.  Computations  demonstrate 
that  these  algorithmic  components  are  also  crucial  to  the  success  of  interface  treatments  created  in  multi- 
domain  and  domain-decomposition  strategies.  For  three-dimensional  computations,  special  metric  relations 
are  employed  to  assure  the  fidelity  of  the  scheme  in  highly  curvilinear  meshes.  A  variety  of  problems, 
including  several  benchmark  computations,  demonstrate  the  success  of  the  overall  computational  strategy. 

Chapter  5  documents  the  critical  extension  of  the  methodology  to  overset  grid  systems.  Overset  grids 
allow  the  efficient  and  accurate  representation  of  wave  propagation  of  the  complex  geometry  often  associated 
with  Air  Force  systems  and  their  components.  A  new  preprocessing  code  BELLERO  has  been  developed  to 
automate  many  of  the  tasks  associated  with  domain  decomposition  for  the  parallel,  high-order  overset-grid 
(HO-OG)  flow  solver  FDL3DI.  The  previous  approach  required  considerable  user  involvement  as  well  as 
manual  modifications  to  the  code  to  set  up  problems  for  processing  using  the  parallel  HO-OG  algorithm. 
Highlighted  capabilities  of  BELLERO  include;  (1)  automatic  generation  of  the  grid  indices  for  the  domain  de¬ 
composition,  taking  into  account  minimum  stencil  requirements  for  the  high-order  algorithm,  (2)  automatic 
generation  of  the  block-level  connectivity  including  periodic  boundary  conditions,  (3)  automated  decom¬ 
position  of  the  grid-level  boundary  conditions,  thus  eliminating  the  need  to  manually  specify  block-level 
boundaries  in  the  code,  and  (4)  calculation  of  high-order  interpolation  coefficients  and  management  of  hole 
points  to  fully  implement  the  HO-OG  approach.  Improvements  have  also  been  made  to  the  FDL3DI  solver 
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itself  in  order  to  further  enhance  the  overall  flexibility  of  the  HO-OG  implementation.  The  new  capability 
is  validated  using  the  benchmark  problems  of  acoustic  electromagnetic  scattering. 

Chapter  6  documents  the  culminating  cleomonstration  of  the  high-order  methodology  on  a  critical  Air 
Force  challenge:  radio  frequency  propagation  through  ionized  gases.  Re-entry  systems  are  engulfed  in  a 
plasma  that  interferes  with  the  communication  signals  to  and  from  the  vehicle.  This  section  includes  details 
of  the  research  required  to  characterize  the  distribution  of  ions  in  the  plasma  sheath  and  the  construction 
of  detailed  grid  system  to  represent  both  the  flight  vehicle  and  the  horn  antennae.  The  high-order,  overset 
method  was  applied  to  the  RAM-C  re-entry  configuration  and  flight  conditions  and  representative  results 
for  the  wave  propagation  response  and  detailed. 

This  report  describes  the  development  and  successful  demonstration  of  high-order  methods  for  wave 
propagation  for  Air  Force  relevant  challenges.  In-house  researchers  in  the  Computational  Sciences  Branch 
of  the  Air  Force  Research  Laboratory  will  exploit  this  new  technology  and  extend  the  work  in  this  area 
in  pursuit  of  other  Air  Force  priorities.  The  next  phase  of  investigation  will  focus  on  the  impact  of  aero- 
optic  degradation  for  directed  energy  weapons,  using  the  tools  and  methods  developed  under  this  concluding 
prioject. 
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Chapter  2 

INTRODUCTION 


An  important  aspect  of  both  civilian  and  military  aircraft  operation  is  the  impact  of  aerodynamically  gen¬ 
erated  sound  on  communities  and  structures.  Examples  of  applications  of  current  interest  include  weapon 
cavity  acoustics,  jet  screech,  sonic  boom,  cabin  noise,  and  sound  generated  by  blade/vortex  interactions. 
In  particular,  the  need  to  meet  more  stringent  community  noise  level  standards  has  resulted  in  increased 
attention  being  given  to  the  relatively  new  field  of  time-domain  computational  aeroacoustics  (CAA).  CAA 
focuses  on  the  accurate  prediction  of  aerodynamic  sound  generated  by  airframe  components  and  propulsion 
systems,  as  well  as  on  its  propagation  and  far-held  characteristics.  Both  aspects  of  the  problem  (i.e.,  sound 
generation  and  propagation)  are  extremely  demanding  from  a  time-domain  computational  standpoint  due 
to  the  large  number  of  grid  points  and  small  time-steps  that  are  typically  required.  Therefore,  for  realistic 
aeroacoustic  simulations  to  become  more  feasible,  higher-order  accurate  and  optimized  numerical  schemes 
are  sought  in  order  to  reduce  the  required  number  of  grid  points  per  wavelength  while  still  ensuring  tolerable 
levels  of  numerically-induced  dissipation  and  dispersion.  These  strict  simulation  requirements  are  shared  by 
other  time-domain  linear  wave  propagation  disciplines  such  as  computational  electromagnetics. 

Recent  reviews  of  computational  aeroacoustics  have  been  given  by  Tam  [2]  and  Wells  and  Renaut  [3]  who 
discuss  various  numerical  schemes  currently  popular  in  CAA.  These  include  among  others  the  dispersion- 
relation-preserving  (DRP)  scheme  of  Tam  and  Webb  [4],  the  method  for  minimization  of  group  velocity 
errors  (MGV)  due  to  Holberg  [5],  the  family  of  high-order  compact  differencing  schemes  of  Lele  [6]  and 
Essentially  Non-Oscillatory  (ENO)  schemes  [7].  The  DRP,  MGV  and  compact  schemes  are  all  centered 
non-dissipative  schemes,  a  property  which  is  desirable  for  linear  wave  propagation.  However,  the  inherent 
lack  of  numerical  dissipation  may  also  result  in  spurious  numerical  oscillations  and  instability  in  practical 
applications  involving  general  geometries,  approximate  boundary  conditions  or  non-linear  features.  In  the 
DRP  approach  for  instance,  artificial  selective  damping  must  be  employed  under  these  conditions  [2].  While 
quite  robust,  standard  upwind  and  upwind-biased  formulations  may  be  undesirable  for  situations  involving 
linear  wave  propagation  due  to  their  excessive  dissipation.  To  overcome  this  difficulty,  higher-order  upwind  [8] 
or  ENO  approaches  must  be  employed.  The  above  spatial  semi-discretizations  are  typically  combined  with 
high-order  explicit  time-integration  methods  such  as  the  multi-stage  Runge-Kutta  procedure.  In  addition 
to  the  spatial  and  temporal  discretizations,  another  critical  aspect  in  CAA  simulations  is  the  accurate 
treatment  of  the  physical  and  computational  boundary  conditions.  A  recent  review  of  radiation,  outflow  and 
wall  boundary  treatments  has  been  provided  in  Ref.  [9]. 

This  work  effort  focused  on  the  development  and  evaluation  of  a  high-order  computational  methodology 
for  aeroacoustic  and  electromagnetic  simulation  on  general  geometries.  There  are  two  primary  components 
to  the  algorithm  chosen  in  the  present  work.  The  first  refers  to  the  differencing  scheme,  for  which  the  choice 
rests  primarily  on  Pade-type  tridiagonal-based  fourth-  and  sixth-order  compact-difference  formulas  [6,  10]. 
As  compared  to  explicit  schemes,  this  approach  incurs  an  increase  in  computational  cost  but  lowers  the 
dispersion  error  and  reduces  the  number  of  points  near  the  boundary  where  special  formulations  are  required. 
As  noted  earlier,  the  non-dissipative  nature  of  centered  schemes  make  them  susceptible  to  the  unrestricted 
growth  of  spurious  perturbations.  The  remedy  employed  here  is  based  on  filtering,  which  comprises  the 
second  and  equally  important  component  of  the  overall  scheme.  The  filters  considered  in  the  present  work 
are  based  on  Pade-type  formulas  requiring  the  solution  of  tridiagonal  systems  of  equations  in  the  general 
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case.  The  expressions  are  taken  from  Ref.  [1],  where  up  to  lOth-order  filters  were  employed  to  stabilize 
finite-volume  electromagnetics  calculations,  and  have  been  successfully  extended  to  solve  the  Navier-Stokes 
equations  [11,  12,  13,  14].  The  filtering  strategies  of  Ref.  [11]  were  also  shown  to  be  promising  for  simple 
aeroacoustic  simulations  on  stretched  Cartesian  meshes  [15]. 

To  accurately  simulate  complex  configurations,  overset-grid  methods[16]  (also  referred  to  as“Chimera” 
methods)  are  widely  used  today.  The  approach  typically  grids  the  entire  computational  domain  using 
multiple  overlapping  structured  component  meshes.  The  use  of  overlapping  component  grids  relieves  many 
of  the  inherent  limitations  associated  with  structured  grids  when  dealing  with  complex  geometries.  Overset 
grid  methods  have  been  shown  to  be  an  effective  way  of  dealing  with  complex  geometries [17,  18],  moving- 
body  problems[19],  or  providing  a  grid-adaptation  capability  [20].  A  summary  of  current  overset-grid  research 
for  a  variety  of  applications  may  be  found  in  Ref.  [21]. 

High-orcler  finite-difference  methods[22,  23]  for  general  curvilinear  grids[24]  using  Pade-type  formulations 
for  the  differencing  and  filtering  operations  [6]  have  also  been  employed  to  perform  numerical  simulations  for 
a  variety  of  problems.  This  approach  has  been  shown  to  provide  highly  accurate  results  for  computations 
involving  turbulent  flows [25,  26],  aeroacoustics[27],  electromagnetics [28,  1]  and  magnetogasdynamics[29]. 
Recently,  this  algorithm  has  been  extended  to  handle  general  overset  grids  in  order  to  provide  additional 
flexibility  to  the  algorithm  for  handling  more  complex  geometries  and  to  provide  localized  grid  refinement.  [30] 
This  high-order  overset-grid  (HO-OG)  approach  has  been  validated  against  fundamental  benchmark  problems 
in  turbulent  flow[31],  aeroacoustics[32]  and  electromagnetics [33]. 

Another  longstanding  problem  for  the  military  is  the  radio-blackout  problem  for  re-entry  vehicles.  At 
high  Mach  numbers,  the  flow  behind  the  shock  can  become  hot  enough  to  ionize  the  flow.  This  has  long 
caused  problems  with  communication  and  has  been  an  active  area  of  research  for  nearly  half  a  century.  The 
problem  has  been  difficult  to  address  because  of  the  complexity  of  the  chemistry  during  re-entry  as  well 
as  the  disparate  length  scales  between  the  wavelength  of  the  signal  and  the  dimensions  of  the  vehicle.  A 
number  of  studies  of  transmission  blackout  were  studied  in  the  1960s  using  rocket  launched  experimental 
probes.  These  tests,  known  as  the  Radio  Attenuation  Measurement  or  RAM-C  tests  attempted  to  measure 
the  electron  density  and  signal  attenuation  at  hypersonic  speeds  [34,  35].  Recently  progress  has  been  made 
in  simulating  these  sorts  of  non-equilibrium  flows  by  Josyula  and  Bailey  [36]. 

For  simulating  wave  propagation  in  dispersive  media,  a  common  method  is  to  use  an  auxiliary  difference 
equation  approach  [37].  For  simulating  wave  propagation  through  weak  plasmas,  a  more  physical  model 
of  the  difference  equation  has  been  used  with  4th-order  FDTD  by  Young  [38]  and  more  recently  a  similar 
formulation  with  3rd-orcler  MUSCLE  by  Shang  [39].  The  current  effort  has  applied  the  technology  developed 
in  simulation  for  acoustics  and  electromagnetics  as  well  as  the  tools  for  overset-grid  technology  to  investigate 
this  problem. 
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Chapter  3 


METHODS  AND  PROCEDURES 


3.1  Governing  Equations 


In  order  to  develop  a  procedure  suitable  for  nonlinear  aeroacoustic  and  plasma  applications  over  complex  con¬ 
figurations,  the  equations  are  cast  in  strong  conservative  form  by  introducing  a  general  curvilinear  coordinate 
transformation  ( x,y,z )  — >  (£,77)0-  In  vector  notation,  these  equations  are: 


d  /U\  <9F  dG  <9H  _  S 

Ft\7 J  +  d£  +  ~d^  +  ac  “  J 


(3.1) 


Here  U  denotes  the  solution  vector;  J  is  the  transformation  Jacobian;  F,G,H  are  the  fluxes  and  S  is  a 
vector  source  term. 


3.2  Numerical  Procedure 

A  finite-difference  approach  is  employed  to  discretize  the  above  equations.  All  discrete  quantities  are  therefore 
assumed  to  be  pointwise  in  nature.  This  choice  is  motivated  by  the  relative  ease  of  formal  extension  to  higher- 
order  accuracy. 


3.2.1  Differencing  Scheme 

For  any  scalar  quantity,  </>,  such  as  a  metric,  flux  component  or  flow  variable,  the  spatial  derivative  (j)'  is 
obtained  in  the  transformed  plane  by  solving  the  tridiagonal  system: 


<x4>'i- 1 


+  4>'i  +  a(t)'i+ i 


4>i+ 2  -  4>i-  2 
b - 4 - 


—  (fri  —  2 

~2 


(3.2) 


where  a,  a  and  6  determine  the  spatial  properties  of  the  algorithm.  The  formula  yields  the  compact  five- 
point,  sixth-order  06,  and  three-point  fourth-order  04  schemes  with  a  =  1/3,  a  =  14/9,  6=1/9  and  a  = 
1/4,  a  =  3/2,  6  =  0  respectively.  Equation  3.2  also  incorporates  the  standard  explicit  fourth-order  E4  ( a 
=  0,  a  =  4/3  and  6  =  -1/3)  and  second-order  E2  (a  =  0,  a  =  1,  6  =  0)  schemes.  At  boundary  points  1,  2, 
IL  —  1  and  IL,  higher-order  one-sided  formulas  are  utilized  which  retain  the  tridiagonal  form  of  the  equation 
set.  These  are  described  in  more  detail  in  Ref.  [40]. 

The  derivatives  of  the  fluxes  are  obtained  by  first  forming  these  fluxes  at  the  nodes  and  subsequently 
differentiating  each  component  with  the  above  formulas. 


3.2.2  Metric  Evaluation 

The  extension  of  high-order  schemes  to  non-trivial  3-D  geometries  demands  that  issues  of  freestream  preser¬ 
vation  and  metric  cancellation  be  carefully  addressed.  These  errors,  which  arise  in  finite-difference  discretiza¬ 
tions  of  governing  equations  written  in  strong-conservation  form,  can  catastrophically  degrade  the  fidelity 
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of  higher-order  approaches.  In  Ref.  [11],  it  was  shown  that  on  highly  distorted  curvilinear  2-D  meshes,  the 
compact  scheme  exhibits  very  small  metric  cancellation  errors  when  the  metrics  are  evaluated  with  the  same 
finite-difference  expressions  as  those  employed  for  the  fluxes.  This  work  also  clearly  showed  that  the  practice 
of  prescribing  analytic  metrics  on  stretched  curvilinear  meshes  can  lead  to  unacceptable  errors  and  therefore 
should  in  general  be  avoided. 

The  situation  is  far  different  in  3-D  however.  As  discussed  in  Refs.  [41]  and  [42]  for  the  second-order 
scheme,  the  previous  straightforward  approach  of  calculating  the  metrics,  although  effective  in  2-D,  fails  to 
provide  metric  cancellation  for  general  3-D  curvilinear  meshes.  To  illustrate  this  point,  consider  the  metric 
relations: 


€x/J  =  yr,ZC  -  VcZr, 

Vx/J  =  vczt  -  vtz c  (3-3) 

Cx/ J  —  y^zr]  yr]Z^ 

Associated  with  this  component  of  the  surface  area  vectors,  there  is  the  corresponding  metric  identity  [42] 

+  {flx/J)ri  +  {Cx/J)c  =  0  (3.4) 

which  must  be  satisfied  numerically  to  ensure  freestream  preservation.  Similar  relations  exist  for  the  other 
two  components  of  the  surface  area  vectors.  Evaluation  of  the  y  and  2  derivatives  in  Eqn.  3.3  using  explicit 
or  compact  centered  schemes  does  not  satisfy  the  identity  of  Eqn.  3.4.  Therefore  grid-induced  errors  may 
appear,  for  instance  in  regions  of  large  grid  variation  or  near  singularities.  Pulliam  and  Steger  [41]  introduced 
a  simple  averaging  procedure  which  guarantees  freestream  preservation  on  general  3-D  curvilinear  meshes. 
This  approach,  which  works  very  well  for  the  2nd-order  scheme,  is  not  readily  extendable  to  higher-order 
formulations.  An  alternate  method  to  enforce  the  metric  identities  consists  in  rewriting  the  metric  expressions 
of  Eqn.  3.3  prior  to  discretization  in  the  equivalent  “conservative”  form[42]: 

£x/J  =  {yrjz)(  ~  i.y(z)r] 

Vx/J  =  (z/c2)?  -  (V£z) C  (3-5) 

Cx/J  =  {yt;z)r,  -  (yr,z)i 

with  similar  relations  for  the  remaining  metric  terms.  When  the  metric  expressions  of  Eqn.  3.5  are  adopted, 
and  the  derivatives  are  evaluated  with  the  same  high-order  formulas  employed  for  the  fluxes,  freestream 
preservation  is  again  recovered  in  general  3-D  curvilinear  grids[12]. 


3.2.3  Time  Integration 

The  equations  are  integrated  in  time  with  the  classical  fourth-order  four-stage  Runge-Kutta  method.  With 
R  denoting  the  residual,  the  governing  equation  is: 


6>U  (d  F  dG  <9H  S\ 

~di  ~  ~~  +  ~d^  +  ~d(  J ) 


(3.6) 


The  classical  four  stage  method  integrates  from  time  t0  (step  n)  to  t0  +  At,  (step  n+1)  through  the  operations: 

k0  =  At  R(  U0)  h  =  At  R(  Ui) 

k2  =  At  R(  U2)  k3  =  At  R(  U3)  (3.7) 

U”+i  =  Vn  +  -(k0  +  2k i  +  2  k2  +  kS) 

6 

where  U0  =  Ura  =  U(x,  y,  z,  to),  Ui  =  U0  +  fco/2,  U2  =  U0  +  fci/2  and  U3  =  U0  +  k2.  The  scheme  is 
implemented  in  the  low  storage  form  described  in  Ref.  [43]  requiring  three  levels  of  storage. 
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3.2.4  Scattered  Field  /  Total  Field  formulation 

For  many  problems  in  (linear)  acoustics  and  electromagnetics,  it  is  desirable  to  be  able  to  allow  the  incident 
waves  to  propagate  from  outside  the  domain.  In  order  to  model  these  waves,  the  solution  is  divided  into  what 
is  known  as  total  field  and  scattered  field  zones(TFZ/SFZ)[44,  45,  46].  In  the  typical  TFZ/SFZ  interface, 
the  incident  field  is  added  or  subtracted  from  the  solution  when  the  zonal  interface  is  crossed.  Consider  the 
following  discrete  field  distribution: 

U,  -  <  ,if  (3.8) 

Uj  =  U*ot,Vj  >  jf  (3.9) 

where  jf  is  the  interface  location  located  between  two  grid  points  and  \Jtot  =  U scat  +  Umc.  So  for  a  simple 
central  difference  located  at  jf  ft  \  would  give: 


(3.10) 


dTU 


total  I 


U tot2[  -  Usc“‘i  -  \Jinc(x  =  Xj  _4) 
Of'  2 _ Of  2 _  Of  ° 

2Ax 


(3.11) 


However,  for  high-order  schemes,  addition  of  this  interface  requires  extensive  modification  to  the  compact 
stencils[47].  Instead,  the  current  parallel  version  of  the  code  utilizes  the  multiple  grid  paradigm  to  accomplish 
the  zonal  interface.  In  this  formulation,  each  grid  is  designated  as  either  a  total  field  grid  (ZONE=l)  or  a 
scattered  field  grid  (ZONE=0).  When  the  data  is  received  from  another  grid,  it  is  modified  by: 


Field  =  Field  +  (ZON  E(current  grid)  —  ZONE(donor  grid))  *  [Incident  Field]  (3.12) 


Figure  3.1  shows  the  flexibility  of  the  zonal  interfaces  using  grid. 

As  the  imposed  incident  field  is  an  analytical  function  in  time,  it  has  been  shown  that  implementation 
of  the  field  in  the  standard  form  for  Runge-Kuttas  can  introduce  errors  that  reduce  the  temporal  accuracy 
of  the  solution[48,  47].  In  standard  form,  the  Runge-Kutta  stages  would  be: 

k0  =  AtR(Vinc(t0),V0)  h  =  At  R(Uinc(t0  +  At/2),  Ur) 

k2  =  AtR(Uinc(t0  +  At/2),U2)  k3  =  AtR{Umc(t0  +  At),U3)  (3.13) 

U”+i  =  Un  +  l(k0  +  +  2  k2  +  kS) 

6 

However,  for  linear  problems  the  errors  can  be  reduced[48,  47]  if  the  incident  field  functions  are  replaced  by 
the  following  analytical  expression  to  match  the  truncation  terms  of  the  individual  steps: 


U*  =  Uinc(t0) 


Umc 


Ufc  =  U»c(to)  +  A‘a|!!(to) 

Uinc(to)  +  f ^T(io)  +  ^^F(to) 


umc  =  U  inc(t0)  +  At 


duin 

dt 


'-(to) 


At 2  9JU!" 
2  dt2 


~(to) 


a F 

4  dt3 


~(to) 


(3.14) 


3.2.5  Interior  filtering  scheme 

The  filtering  procedure  forms  an  important  component  of  the  present  centered  algorithm  since  its  function 
is  to  suppress  numerical  instabilities  arising  from  unresolved  scales,  mesh  non-uniformities  and  boundary 
conditions.  If  a  typical  component  of  the  solution  vector  is  denoted  by  (f>,  filtered  values  tp  satisfy, 

Otf4>i- 1  +  <f>i  +  CZf<t>i- 1-1  =  Q—  (4>i+n  +  4>i-n)  (3.15) 
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Table  3.1:  Coefficients  for  filter  formula  at  interior  points  [1].  ay  is  a  free  parameter  in  the  range  0  <  |ay|  < 
0.5. 


Scheme 

a0 

ai 

a2 

a3 

a  5 

OA 

F2 

h  +  af 

k  +  af 

0 

0 

0 

0 

2 

F4 

5  ,  3a/ 

8  '  4 

\+af 

-1  1  af 

8  _r  4 

0 

0 

0 

4 

F6 

11  ,  5a/ 

16  '  8 

15  ,  17a/ 

32  '  16 

-3  ,  3a/ 

16  '  8 

1 

32  16 

0 

0 

6 

F8 

93+700/ 

128 

7+18a/ 

16 

—  7+14a  / 

32 

1  af 

16  8 

-1  1  <*f 

128  '  64 

0 

8 

F10 

193+126a/ 

105+302a  / 

15(-l+2a/) 

45(1  — 2a/) 

5(  —  l+2a/) 

1— 2a/ 

10 

256 

256 

64 

512 

256 

512 

For  multi-dimensional  problems,  the  filter  is  applied  sequentially  in  each  of  the  three  directions.  Equation  3.15 
is  based  on  templates  proposed  in  Refs.  [6]  and  [49]  and  with  proper  choice  of  coefficients,  provides  a  21Vth- 
orcler  formula  on  a  2N  + 1  point  stencil.  The  N  + 1  coefficients,  a 0,  a\, . . .  a jv,  are  derived  in  terms  of  ay  with 
Taylor-  and  Fourier-series  analyses  and  are  listed  in  Table  3.1.  On  uniform  meshes,  the  resulting  filters  are 
non-dispersive,  do  not  amplify  any  waves,  preserve  constant  functions  and  completely  eliminate  the  odd-even 
mode.  Since  ay  is  a  free  parameter,  an  explicit  filter  i.e.,  not  requiring  the  solution  of  a  tridiagonal  can 
be  easily  extracted  by  setting  ay  =  0.  In  this  case,  the  formulas  of  Ref.  [50]  are  identically  recovered.  The 
primary  constraint  on  ay  is  that  it  must  satisfy  the  inequality  —0.5  <  ay  <  0.5.  In  this  range,  higher  values 
of  ay  correspond  to  a  less  dissipative  filter,  and  at  ay  =  0.5,  there  is  no  filtering  effect.  Detailed  spectral 
responses  of  these  filters  may  be  found  in  Refs.  [1]  and  [40].  For  the  results  below,  the  solution  is  filtered 
once  (in  each  spatial  direction)  after  the  final  stage  of  the  explicit  RK4  scheme.  The  interior  filtering  formula 
utilized  is  denoted  by  appending  its  designation  to  that  of  the  scheme.  For  example,  C6F10  designates  the 
sixth-order  compact  scheme  combined  with  the  tenth-order  filter. 


3.2.6  Near-boundary  filter  formulations 

The  relatively  large  stencil  of  high-order  filters  requires  special  formulations  at  several  points  near  the 
boundaries.  For  instance,  the  lOth-orcler  interior  filter  requires  an  11  point  stencil  and  thus  cannot  be 
applied  at  the  “near-boundary”  points  1, . . . ,  5  and  correspondingly  at  IL  —  4, . . . ,  IL  where  it  protrudes  the 
boundary.  The  values  at  points  1  and  IL  are  specified  explicitly  through  the  boundary  conditions  and  are 
not  filtered.  At  the  remaining  near-boundary  points,  the  approach  followed  in  Ref.  [11]  was  to  reduce  the 
stencil-size  by  applying  lower-order  centered  ( LOC)  formulas.  Thus,  for  example,  at  points  2  and  IL  —  1,  a 
2nd-order  filter  is  applied,  at  3  and  IL  —  2,  a  4th-order  filter  is  applied  and  so  on.  However,  as  the  mesh  is 
coarsened  near  the  boundaries,  the  error  induced  by  this  low-order  central  or  LOC  technique  may  eventually 
become  unacceptable  and  adversely  affect  the  global  solution  accuracy.  In  such  instances,  the  lower-order 
filter  at  points  2,3,  IL  —  2,  IX  —  1  may  be  optimized  by  specifying  the  filter  control  parameter  ay  to  be  as 
close  as  possible  to  0.5  as  stability  allows.  This  optimized  low-order  central  or  OLOC  approach  was  shown 
to  be  effective  for  the  solution  of  viscous  flows  in  Ref.  [11]  and  clearly  illustrates  an  advantage  of  the  present 
implicit  filters  over  standard  explicit  formulations  which  lack  control  of  the  filter  spectral  response. 

A  more  general  alternative  to  the  OLOC  approach,  is  the  use  of  the  higher-order  one-sided  filter  formulas 
developed  in  Ref.  [12].  At  a  near-boundary  point,  i,  a  filter  formula  is  given  by 


Qftfti—l  4“  4“  OLf(j>i-\-\  —  Yjn=^Cln,i4>n 

i  e  {2,.. .,5} 

1  4“  4*i  4”  1-1  ^-2n=0®'IL  —  n,i4*I L— m 

i  e  {IL  —  4, . . . ,  IX  —  1} 


(3.16) 


This  choice  retains  the  tridiagonal  form  of  the  filter,  and  ay  remains  as  the  only  free  parameter. 
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Table  3.2:  Coefficient  for  6i/i-order  boundary  filter  formula  at  point  3 


OA 

an, 3 

02,3 

03,3 

014,3 

05,3 

06,3 

07,3 

6 

1  1  af 

64  '  32 

3  |  13af 

32  '  16 

49  ,  15q:/ 

64  '  32 

5  ,  3af 

16  '  8 

15  , 

64  '  32 

3  3af 

32  16 

1  1  af 

64  '  32 

Table  3.3:  Coefficients  for  boundary  filter  formulas  at  point  2 

OA 

an,  2 

o  2,2 

0-3,2 

04,2 

05,2 

06,2 

07,2 

4 

1  ,  (Cif 

16  '  8 

3  _i_  af 

4  '  2 

3  , 

8  '  4 

1  1  af 

4  '  2 

1  af 

16  8 

0 

0 

6 

1  ,  31af 

64  '  32 

29  i  3af 

32 r  16 

15  ,  17af 

64  '  32 

5  ,  5a/ 

16  '  8 

15  15af 

64  32 

3  i  3af 
32  '  16 

1  af 

64  32 

Figure  3.2:  Spectral  function  of  some  near-boundary  filter  formulas  at  point  2 


Tables  3.2  and  3.3  list  coefficients  for  the  fourth-  and  sixth-order  left-boundary  filter  formulas  employed  in 
the  present  computations  at  points  2  and  3.  The  right-boundary  formulas  are  obtained  by  noting  ajL-n,i  = 
Un+i.iL-i+i  for  i  e  {IL  —  4, . . . ,  IL  —  1}.  An  extensive  listing  of  the  boundary  filter  coefficients  is  provided 
in  Refs.  [12,  40].  The  real  component  of  the  spectral  function  SF  of  some  boundary  formulations  at  point  2 
are  plotted  in  Fig.  3.2.  The  higher-order  approaches  are  observed  to  be  superior  low-pass  filters  in  the  range 
of  interest,  w  <  1.2  corresponding  to  roughly  5  or  more  points  per  wave.  The  2nd-order  filter  is  symmetric 
at  point  2  and  its  real  part  is  always  less  than  one,  i.e.,  no  wave  is  amplified.  The  asymmetric  higher-order 
formulas  do  not  satisfy  this  desirable  Real(SF)  <  1  constraint.  The  range  of  amplification  increases  with 
order  of  accuracy  (one-sidedness  of  the  formula)  and  for  a  given  order  of  accuracy,  the  implicit  (a/  ^  0) 
filter  is  seen  to  be  better  than  the  explicit  variant  (a/  =  0)  in  this  respect.  The  near-boundary  filter  spectral 
response  may  be  optimized  even  further  by  increasing  a.f  thereby  reducing  the  undesirable  amplification 
behavior.  This  is  demonstrated  in  Fig.  3.2  for  the  6th-order  filter  with  a/  =  0.49.  This  advantage  of  the 
implicit  filter  will  be  illustrated  in  the  scattering  results  presented  below. 
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Chapter  4 

COMPUTATIONAL 

AEROACOUSTICS 

4.1  Introduction 

The  numerical  procedures  given  in  the  previous  chapter  have  been  applied  to  the  simulation  of  aeroacoustic 
phenomena.  The  equation  set  and  boundary  conditions  are  presented  here  followed  by  a  number  of  test  cases 
to  validate  the  methodology. 


4.2  Governing  Equations 

To  accurately  simulate  nonlinear  aeroacoustic  applications,  the  full  (non-linear)  Euler  equations  are  selected 
in  general  curvilinear  coordinates.  As  mentioned  in  the  previous  chapter,  the  form  of  the  equations  may  be 
written  as: 

d  (u\  dF  dG  dH  S 

at  7  +aTfai+a<=J  (41) 


Here  U  =  {p,  pu ,  pv,  pw,  pE}  denotes  the  solution  vector,  J  is  the  transformation  Jacobian,  and  F.  G  and  H 
are  the  inviscid  fluxes: 

pU 

puU  +  £xp 

pvU  +  CyP  I  (4.2) 

pwU  +  £zp 
(. pE  +  p)U 


where 


F  = 


G  = 


pV 

puV  +  TJxP 
pvV  +  rjyp 
pwV  +  r/zp 
(pE  +  p)V 

pW 

puW  +  CxP 
pvW  +  C yP 
pwW  +  ( zp 
(, pE  +  p)W 


U  =  £,XU  +  CyV  +  £ZW 

V  =  PxU  +  rjyV  +  rjzw 


(4.3) 


(4.4) 
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(4.5) 

(4.6) 


W  =  (XU  +  CyV  +  c zW 


E  = 


T 


+  \{u2  +  v2  +  w2). 


(4.7) 

(4.8) 


(7-1  )Ml  '  2' 

In  the  expressions  above,  u,  v,  w  are  the  Cartesian  velocity  components,  p  the  density,  p  the  pressure,  and  T 
the  temperature.  The  perfect  gas  relationship  p  =  pRT  is  also  assumed. 

The  vector  source  term  S  =  {0, 0,  0,  0,  S$}  is  included  in  order  to  account  for  acoustic  sources  which,  for 
the  present  results,  are  prescribed  to  be  of  the  form: 


^  /  X  lno  (x-xc)2  +  (y-yc)2 

£5  (x,y,t)  =  e  1  &  sin  (u>t)  f  {t) 


f  ( t )  =  min  1.0, 


(4.9) 


where  xc  and  yc  denote  the  center  of  the  source,  u>  is  the  frequency  and  f(t)  is  the  function  used  to  ramp 
the  onset  of  the  source  at  the  beginning  of  the  computations.  The  parameters  in  Equation  4.9,  including  b 
and  ta,  are  specified  below  in  the  description  of  the  results. 


4.2.1  Boundary  Conditions 

In  the  present  acoustic  simulations,  several  types  of  boundaries  occur  where  wall,  symmetry  or  farfield 
(outflow)  conditions  need  to  be  specified.  At  walls,  standard  inviscid  conditions  are  employed,  in  which 
the  normal  component  of  velocity  is  set  to  zero  and  the  pressure,  density  and  tangential  components  of 
velocity  are  extrapolated.  All  Neumann-type  conditions,  including  those  required  to  enforce  symmetry,  are 
approximated  to  higher-order  accuracy  with  one-sided  third-  or  fourth-order  formulas. 

Although  several  sophisticated  treatments  are  possible  for  absorbing  or  non-reflecting  farfield  condi¬ 
tions  [9] ,  the  present  work  employs  a  simple  yet  robust  technique  developed  by  exploiting  the  characteristics 
of  the  low-pass  filter.  As  discussed  in  more  detail  in  Section  4.1.3,  rapid  stretching  in  regions  outside  the 
domain  of  interest,  promotes  damping  of  all  perturbations  via  transfer  of  energy  to  reflected  odd-even  modes 
which  are  in  turn  annihiliated  by  the  baseline  filter.  Although  the  present  approach  is  empirical  in  nature, 
it  offers  the  potential  of  extension  to  non-linear  and  multi-dimensional  situations  where  techniques  derived 
through  asymptotic  or  linear  analyses  may  be  ineffective. 

4.2.2  Interface  treatment  in  multi-domain  calculations 

A  major  issue  in  multi-domain  computations  is  the  accurate  treatment  of  interfaces  between  domains.  In  the 
present  work,  communication  between  adjacent  meshes  is  conducted  through  finite-size  overlaps,  as  depicted 
in  the  schematic  of  Fig.  4.1  for  a  five-point  vertical  overlap.  At  every  time-step,  the  solution  is  advanced 
independently  in  each  domain  with  individual  interior  and  boundary  formulas  in  the  same  manner  as  in 
single-domain  computations.  Data  is  exchanged  between  adjacent  domains  at  the  end  of  each  stage  of  RK4, 
as  well  as  after  each  application  of  the  filter.  Each  vertical  line  is  denoted  by  its  i-index.  The  values  at  points 
1  and  2  of  Mesh  2  are  set  to  be  identically  equal  to  the  corresponding  updated  values  at  points  IL  —  4  and 
IL  —  3  of  Mesh  1.  Similarly,  reciprocal  information  is  transferred  through  points  4  and  5  of  Mesh  2  which 
“donate”  values  to  points  IL  —  1  and  IL  of  Mesh  1.  This  exchange  of  information  is  shown  with  arrows 
at  each  point  in  the  schematic.  Larger  and  smaller  overlap  regions  have  been  examined  for  fluid  dynamic 
applications  in  Ref.  [12]. 


4.3  Results 


4.3.1  Impact  of  Filtering  on  Stretched  Meshes 

In  practical  scattering  simulations,  grid  stretching  is  usually  employed  in  order  to  reduce  required  com¬ 
putational  resources  as  well  as  to  permit  the  use  of  approximate  farfield  boundary  conditions.  Therefore, 
the  performance  of  high-order  schemes  on  general  stretched  meshes  must  be  carefully  examined  for  the 
methodology  to  be  applicable  to  complex  configurations. 
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Mesh  1 


O  Mesh  2 


Arrows  show 
information  transfer 


Figure  4.1: 


Schematic  of  vertical  5-point  mesh  overlap  for  multi-domain  applications 


One-Dimensional  Pulse 

An  enlightening  analysis  of  the  behavior  of  a  smooth  solution  as  it  passes  through  a  sudden  mesh  coarsening 
has  been  presented  by  Vichnevetsky  [51]  for  the  1-D  advection  equation  semi-discretized  with  the  standard 
second-order  E2  centered  scheme.  This  analysis  indicates  that  although  total  energy  is  preserved,  at  the  grid 
coarsening  interface  a  significant  portion  of  the  energy  is  deposited  on  a  reflected  solution  composed  primarily 
of  odd-even  modes  and  modulated  by  a  smooth  envelope.  This  reflected  energy  propagates  upstream  (i.e. 
with  negative  group  velocity)  and  in  most  circumstances  if  left  unchecked  it  has  the  potential  of  ultimately 
contaminating  the  genuine  solution. 

To  highlight  the  effect  of  grid  stretching,  we  solve  the  1-D  Euler  equations  for  the  propagation  of  a  density 
disturbance  through  the  sudden  grid  coarsening  shown  in  Fig.  4.2.  The  initial  conditions  are  given  by 


u  =  l,v  =  0,p  =  poo 

p  =  1  +  ee-^2-*^ 


(4.1) 


where  e  =  0.1.  It  should  be  noted  that  no  special  one-sided  treatments  are  invoked  at  the  grid  interface, 
but  rather  the  grid  jump  is  handled  through  the  metrics  in  the  general  coordinate  transformation.  Results 
with  the  E2  scheme  are  shown  in  Fig.  4.3.  The  figure  indicates  that  in  accordance  with  Ref.  [51],  when  the 
pulse  crosses  the  interface,  a  reflected  solution  emerges  which  is  characterized  by  spurious  high-frequency 
modes.  Similar  behavior  was  found  for  the  E4,  C4  and  C6  schemes  (not  shown).  At  a  given  time  instant,  the 
upstream  penetration  of  the  reflected  solution  increased  with  stencil  size  and  implicitness  of  the  operator. 

Since  the  reflections  due  to  grid  stretching  are  characterized  by  high-frequency  modes,  they  can  be  easily 
removed  from  the  solution  by  the  high-order  low-pass  filter.  This  is  shown  in  Fig.  4.4(a)  for  the  C6F8 
scheme,  where  spurious  reflections  at  the  grid-coarsening  interface  are  significantly  diminished  by  the  use  of 
the  8th-order  filter  with  a/  =  0.45. 

These  results  suggest  that  the  combination  of  significant  grid  stretching  with  a  discriminating  low-pass 
filter  may  be  used  as  an  alternative  procedure  for  outflow  boundary  treatment.  This  approach  is  investigated 
further  in  Section  4.1.3. 
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Ax1=0.2 


Ax2=2 


Figure  4.2:  Mesh  with  abrupt  change  in  spacing 


Figure  4.3:  Propagation  of  1-D  pulse  through  sudden  mesh  coarsening  with  E2  scheme  and  without  filter. 


Figure  4.4:  Propagation  of  1-D  pulse  through  sudden  mesh  coarsening  with  C6F8  scheme. 
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Two-Dimensional  Pulse 


Next,  we  consider  the  propagation  of  a  2-D  pressure  pulse  in  the  presence  of  grid  stretching.  The  mesh, 
shown  in  Fig.  4.5(a),  contains  a  section  of  uniform  spacing  (Ax  =  Ay  =  0.05)  which  extends  over  the  region 
—3  <  {x,  y}  <  3.  Beyond  this  zone,  the  mesh  is  stretched  with  an  exponential  function  and  maximum 
stretching  factor  of  1.13.  For  the  purpose  of  comparison,  a  large  uniform  grid  was  also  considered.  The 
solutions  were  initialized  at  t  =  0  by  prescribing  a  pressure  disturbance  of  the  form 


P  =  Poo  +  ee 


b2 


(4.2) 


where  e  =  0.1,  b  =  0.2  and  xc  =  yc  =  0.  The  propagation  of  this  pressure  pulse  was  computed  with 
the  C6  scheme  and  At  =  0.005  until  t  =  4.5.  By  this  time,  the  pressure  disturbance  had  moved  into  the 
stretched  mesh  region  but  was  still  sufficiently  far  from  the  computational  farfield  boundaries  in  order  to 
diminish  effects  of  boundary  conditions.  The  pressure  contours  are  displayed  in  Fig.  4.5(b)  and  (c)  for 
both  the  unfiltered  and  filtered  ( F10 ,  a/  =  0.49)  solutions.  The  results  obtained  without  a  filter  show  the 
appearance  of  high-frequency  modes  which  are  generated  as  the  pulse  moves  into  the  stretched  mesh  region. 
As  in  the  previous  1-D  case,  this  spurious  reflected  energy  propagates  in  a  direction  opposite  to  the  travel 
of  the  physical  disturbance,  towards  the  center  of  the  domain.  When  the  high-order  filter  is  employed, 
these  unwanted  reflections  are  completely  eliminated,  as  seen  in  Fig.  4.5(c).  The  corresponding  computed 
pressure  along  the  center  of  the  pulse  ( y  =  0)  is  given  in  Fig.  4.6.  In  the  region  of  interest,  —3  <  x  <  3,  the 
filtered  results  for  the  stretched  grid  are  essentially  identical  to  those  obtained  without  a  filter  on  a  large 
domain  with  constant  spacing  (Ax  =  Ay  =  0.05).  By  contrast,  the  unfiltered  solution  on  the  stretched  grid 
clearly  displays  the  existence  of  packets  of  reflected  energy  composed  mostly  of  odd-even  modes.  These  2-D 
results  clearly  demonstrate  the  effectiveness  of  the  high-orcler  low-pass  filter  in  controlling  spurious  reflections 
without  degrading  the  fidelity  of  the  solution  in  the  region  of  interest.  In  the  absence  of  a  high-frequency 
cutoff  mechanism  (designed  to  completely  annihiliate  odd-even  modes),  a  large  number  of  grid  points  would 
be  required  [15]  if  a  constant  mesh  spacing  had  to  be  carried  all  the  way  to  the  far  field  boundaries. 


Use  of  Mesh  Stretching  and  Filtering  as  an  Alternative  Outflow  Boundary  Treatment 

As  shown  earlier  in  Section  4.1.1,  by  employing  a  large  rate  of  stretching,  a  significant  amount  of  energy  can 
be  reflected  at  the  grid  coarsening  interface.  Provided  this  reflected  energy  is  deposited  into  high-frequency 
modes  (in  the  fine  mesh  region),  it  can  then  be  subsequently  eliminated  by  the  baseline  high-order  low-pass 
filter  without  contaminating  the  genuine  solution.  Furthermore,  since  the  mesh  is  stretched  rapidly  in  the 
buffer  zone,  the  transmitted  energy  is  also  quickly  dissipated  by  the  high-order  filter  (as  the  transmitted 
solution  features  are  represented  by  a  diminishing  number  of  points  per  wave).  This  proposed  method 
eliminates  the  need  for  more  sophisticated  (and  perhaps  more  restrictive)  boundary  conditions  at  the  expense 
of  extending  the  computational  domain.  Although  the  proposed  approach  has  a  mathematical  foundation 
(at  least  based  on  the  1-D  analysis  of  Ref.  [51]),  its  implementation  is  highly  empirical,  and  therefore  its 
utility  must  be  evaluated  in  the  context  of  practical  applications. 

As  a  severe  test  case,  consider  the  propagation  of  acoustic  waves  in  the  grid  of  Fig.  4.7(a).  This  mesh 
is  uniform  (Ax  =  Ay  =  0.05)  in  the  center  of  the  computational  domain  (—3  <  x,y  <  3).  Outside  of  this 
resolved  region,  Ax  and  Ay  are  increased  abruptly  by  a  factor  of  10.  Both  the  case  of  an  acoustic  source 
and  a  transient  pulse  were  computed  with  C6F8  (aj  =  0.45)  in  order  to  examine  the  robustness  of  the 
numerical  approach.  The  source  was  specified  by  Equation  4.9  with  xc  =  yc  =  0,lo  =  5ir,  b  =  0.2  and  ta  =  4. 
A  snapshot  of  the  pressure  is  displayed  in  Fig.  4.7(b).  It  is  apparent  that  the  acoustic  energy  reflected  at 
the  grid-coarsening  interface  is  almost  completely  annihiliatecl  by  the  high-orcler  filter.  The  corresponding 
plot  of  the  instantaneous  pressure  along  the  diagonal  (x  =  y )  is  shown  in  Fig.  4.7(c)  and  indicates  that  the 
transmitted  energy  diffuses  rapidly  in  the  coarse  mesh  region.  For  this  square  grid,  the  circular  waves  cross 
the  grid  interface  at  a  varying  angle  of  incidence  without  apparent  anisotropies  being  introduced. 

For  a  more  quantitative  test,  the  case  of  a  transient  pulse  given  by  Equation  4.2  (with  e  =  0.02,  b  =  0.2) 
was  also  considered.  Figure  4.8(a)  shows  the  history  of  non-dimensional  pressure  at  the  point  x  =  2,y  =  2. 
The  solution  for  the  suddenly  stretched  mesh  is  observed  to  be  in  excellent  agreement  with  results  obtained 
on  a  large  domain  with  uniform  mesh  spacing,  as  demonstrated  by  the  small  difference  between  the  two 
solutions  (Fig.  4.8(b)). 
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Figure  4.5:  Propagation  of  2-D  pressure  pulse  on  smoothly-stretched  Cartesian  mesh:  (a)  grid,  (b)  unfiltered 
C6  and  (c)  filtered  C6F10  results 
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Figure  4.6:  Effect  of  high-order  filter  on  computed  pressure  along  y  =  0  at  t  =  4.5  for  the  2-D  acoustic  pulse 
propagation  of  Fig.  4.5. 


(a) 


Figure  4.7:  Computation  of  2-D  acoustic  source  on  mesh  with  abrupt  stretching,  (a)  Grid  (b)  Pressure 
contours  (c)  Pressure  along  diagonal 
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Figure  4.8:  Pressure  at  x  =  2,  y  =  2  for  2-D  pulse  propagation  on  mesh  with  abrupt  stretching,  a)  Pressure 
history,  b)  discrepancy  between  stretched  and  uniform  grid  solutions 


Grid 

Dimensions 

max 

D 

^@min 

^^max 

G1 

156  x  175 

0.04 

0.5 

1.25 

G2 

252  x  175 

0.04 

0.5 

0.75 

Table  4.1:  Grid  parameters  for  pulse  scattering  problem  in  region  of  interest  (vjD  <  7) 


The  effectiveness  of  the  proposed  outflow  treatment  was  also  evaluated  for  the  convecting  vortex  case 
previously  considered  in  Ref.  [11].  The  computational  domain  of  interest  (—6  <  x,y  <  6)  is  discretized 
using  a  uniform  mesh  spacing  (Ax  =  Ay  =  0.2),  as  shown  in  Fig.  4.9a.  Outside  of  this  region,  the  grid  is 
expanded  employing  a  constant  stretching  factor  of  1.7.  The  vortex  is  convected  along  the  diagonal  of  the 
mesh  in  order  to  provide  a  more  general  test  of  the  outflow  approach.  Contours  of  the  perturbation  velocity 
magnitude  at  various  instants  are  shown  in  Fig.  4.9b-d.  The  vortical  disturbance  is  observed  to  pass  without 
significant  distortion  through  the  corner  of  the  grid-coarsening  interface.  A  comparison  of  the  computed  and 
exact  perturbation  velocity  magnitudes  along  the  diagonal  x  =  y  is  displayed  in  Fig.  4.10  at  t  =  2  and  8. 
The  excellent  agreement  between  the  computed  and  exact  solution  is  clearly  retained  as  the  vortex  crosses 
the  interface. 

4.3.2  Scattering  of  Acoustic  Pulse  From  Cylinder 

In  order  to  validate  the  present  approach  for  curvilinear  geometries,  we  select  as  a  test  case  the  benchmark 
problem  denoted  as  Category  I,  Problem  2  in  the  2nd  CAA  Workshop  of  Ref.  [52].  This  configuration 
(Fig.  4.11)  describes  the  scattering  off  a  circular  cylinder  of  a  prescribed  initial  pressure  pulse.  The  pulse  is 
given  by  Equation  4.2  with  xc  =  4,  yc  =  0,  e  =  0.01,  b  =  0.2. 

Along  the  cylinder  surface,  the  simple  wall  boundary  conditions  noted  in  Section  3.6  are  employed.  Since 
the  configuration  is  symmetric,  only  the  upper  half  of  the  domain  is  considered,  and  symmetry  conditions 
are  invoked  along  6  =  0°,  180°.  As  indicated  in  Fig.  4.11,  the  grid  is  stretched  for  r/D  >  7  using  a  constant 
stretching  factor  of  1.1.  Since  in  this  coarse-mesh  region  the  outgoing  pulse  is  dissipated  by  the  baseline 
filter,  a  simple  extrapolation  condition  is  suitable  along  the  farfield  boundary.  Two  levels  of  spatial  resolution 
were  used,  and  details  of  the  computational  grids  are  provided  in  Table  4.1  .  All  cases  were  advanced  in 
time  with  a  non-dimensional  At  of  0.004. 

Pressure  contours  depicting  pulse  propagation  and  reflection  off  the  cylinder  surface  are  shown  in  Fig.  4.11 
at  several  instants  in  time  for  the  C6F10  scheme  on  the  finest  mesh  (G2).  The  history  of  pressure  at  selected 
points  and  for  several  computations  is  presented  in  Figs.  4.12  through  4.14.  The  points  denoted  as  ‘A’  and 
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Figure  4.9:  Propagation  of  a  vortical  disturbance  at  45°  to  the  horizontal,  a)  Grid.  Contours  of  perturbation 
velocity  magnitude  at  b)  t=2,  c)  t=6  and  d)  t=8 


Figure  4.10:  Perturbation  velocity  magnitude  along  diagonal  for  convecting  vortical  disturbance.  Symbols 
and  lines  denote  computed  and  exact  solutions  respectively 
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Figure  4.11:  Pressure  contours  at  various  instants  for  acoustic  pulse  scattered  by  a  circular  cylinder 


P’ 


Figure  4.12:  History  of  pressure  at  point  ‘A’  on  Grid  G2  for  scattered  pulse 
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p’ 


Figure  4.13:  History  of  pressure  at  point  ‘C’  on  Grid  G2  for  scattered  pulse 


Figure  4.14:  History  of  pressure  at  point  ‘A’  on  Grid  G1  for  scattered  pulse 
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Grid 

Dimensions 

A 

A 

^@max 

T— 5 

At  max 

(t  A0maa; ) 

Gl 

252  x  175 

6.25 

3.82 

0.75 

G2 

361  x  321 

8.00 

5.73 

0.50 

G3 

481  x  491 

12.5 

7.64 

0.375 

Table  4.2:  Grid  parameters  for  source  scattering  problem,  A=  spatial  wavelength 


Figure  4.15:  Pressure  history  at  a  point  for  periodic  acoustic  source  problem 


‘C’  are  located  at  (r/D  =  5, 9  =  90°)  and  ( r/D  =  5, 9  =  180°)  respectively.  Results  on  Grid  G2  with  C4F8 
(af  =  0.4)  are  observed  to  be  in  excellent  agreement  with  the  exact  solution  [52].  This  is  also  the  case  for 
the  C6F10  algorithm  (not  shown).  For  the  purpose  of  comparison,  computations  were  also  performed  with 
a  standard  third-order  accurate  upwind-biased  scheme.  As  Figs.  4.12  and  4.13  indicate,  the  upwind  results 
exhibit  excessive  dissspation  even  on  this  finer  mesh.  On  the  coarser  mesh  Gl,  (Fig.  4.14),  the  C6F1U  and 
C4F8  with  implicit  filter  (af  =  0.4)  are  also  observed  to  be  in  good  agreement  with  the  exact  solution. 
However,  on  this  mesh,  the  C4  scheme  combined  with  an  8th-order  explicit  filter  (i.e.  0/  =  0.0)  displays 
appreciable  error.  This  highlights  the  improved  accuracy  of  the  implicit  filter  formulation. 

4.3.3  Scattering  of  Periodic  Acoustic  Source 

This  case  which  corresponds  to  Category  I,  Problem  1  in  Ref.  [52],  considers  the  scattering  from  a  circular 
cylinder  of  a  periodic  acoustic  source.  An  acoustic  source  with  finite  spatial  support  is  given  by  Equation4.9 
with  xc  =  4,  yc  =  0,  u>  =  87t,  b  =  0.2,  t,0  =  4.  This  case  constitutes  a  more  stringent  test  of  the  algorithm 
and  boundary  conditions  since  an  asymptotic  periodic  solution  must  be  attained,  and  long-term  numerical 
stability  demonstrated.  The  boundary  conditions  were  of  the  same  type  as  those  noted  in  the  previous 
section.  Three  computational  grids  were  employed  and  their  details  are  summarized  in  Table  4.2.  For  all 
grids,  the  radial  spacing  is  stretched  over  the  range  10  <  r/D  <  20  using  a  constant  stretching  factor  of 
approximately  1.2.  A  time  step  At  =  0.004  was  specified  for  all  cases,  and  the  solutions  were  typically  run  for 
approximately  20  —  30  characteristic  times  to  guarantee  a  periodic  state  and  to  ensure  long  term  numerical 
stability. 

A  sample  convergence  history  of  the  pressure  at  r  =  3,6  =  135°  is  shown  in  Fig.  4.15.  After  the  initial 
transient  generated  by  the  source  onset  leaves  the  domain,  a  periodic  solution  is  clearly  obtained.  This 
suggests  that  acoustic  energy  is  not  trapped  by  the  grid  stretching  or  the  approximate  outflow  boundary 
treatment.  A  representative  instantaneous  interference  pattern  generated  by  the  incident  and  scattered 
waves  is  shown  in  Fig.  4.16.  It  is  observed  that  in  the  far  field,  the  rapid  mesh  stretching  combined  with 
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Figure  4.16:  Interference  pattern  due  to  scattering  by  a  cylinder  of  periodic  acoustic  source 
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Figure  4.17:  Computed  directivity  of  radiated  sound  at  r/D  =  5  with  C6F10 


Figure  4.18:  Computed  directivity  of  radiated  sound  at  r/D  =  10  with  C6F10 

the  high-order  filter  provide  an  effective  mechanism  to  absorb  the  outgoing  acoustic  radiation  with  minimal 
reflection.  Figures.  4.17  and  4.18  display  the  computed  directivity  of  the  radiated  sound  D{9)  =  r  <  p'2  > 
at  r/D  =  5  and  r/D  =  10  respectively.  The  directivity  computed  with  C6F10  on  the  finest  mesh  (Grid 
G3 )  is  in  excellent  agreement  with  the  exact  solution  at  both  radial  locations  considered.  The  results  on  the 
medium  mesh  (Grid  G2)  are  reasonable  at  r/D  =  5  but  begin  to  exhibit  errors  at  r/D  =  10,  particularly 
for  0  <  100°.  On  the  coarsest  mesh  (Grid  Gl),  the  computed  D{9)  shows  already  significant  departure  at 
r/D  =  5.  It  should  be  noted  that  the  radial  spacing  in  Grids  Gl ,  G2  and  G3  (see  Table  4.2)  correspond 
to  only  6.25,  8  and  12.5  points  per  wave  (PPW)  respectively.  However,  in  terms  of  the  azimuthal  spacing, 
at  r/D  =  5  the  corresponding  densities  are  only  approximately  4,  6  and  8  PPW.  It  then  appears  that  in 
order  to  maintain  an  accurate  solution  the  minimum  spatial  resolution  requirement  lies  in  the  range  of  4  to  6 
PPW.  For  subsequent  comparisons,  and  in  order  to  limit  use  of  computational  resources,  we  focus  on  Grid 
G2  and  on  the  directivity  at  r/D  =  5  where  the  G6F10  method  provides  good  results. 

A  comparison  of  various  schemes  is  presented  in  Fig.  4.19  in  terms  of  the  computed  directivity.  Limited 
degradation  in  accuracy  is  encountered  when  going  from  the  baseline  G6F10  to  the  C/F8  scheme.  By  contrast, 
the  explicit  fourth-order  method  combined  with  an  explicit  eighth-order  filter  (E/F8)  provides  rather  poor 
results  for  the  sound  pressure  level  at  this  location.  Therefore,  a  significant  improvement  is  achieved  when 
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Figure  4.19:  Effect  of  scheme  on  computed  directivity  at  r/D  =  5 


proceeding  from  the  explicit  to  the  implicit  (compact)  spatial  discretization  and  filtering  techniques  while 
simultaneously  reducing  the  scheme  stencil  size. 

The  effect  of  the  filter  formulation  near  the  scatterer  on  the  solution  accuracy  for  the  C6F10  procedure  is 
examined  next.  Figure  4.20  shows  the  computed  sound  directivity  at  r/D  =  5  for  three  different  boundary 
filter  approaches.  The  high-order  boundary  filter  strategy  (denoted  as  HO)  with  filter  orders  ‘6, 6, 6, 8’  at  points 
‘2, 3, 4, 5’  respectively  is  in  good  agreement  with  the  exact  solution.  With  this  technique,  the  filter  control 
parameter  a/  =  0.45  remains  unchanged  throughout  the  domain.  For  the  low-order-centered  method  (LOC), 
ctf  =  0.45  is  again  fixed,  but  the  filter  order  is  degraded  to  ‘2, 4, 6, 8’  at  points  ‘2, 3, 4, 5’.  For  the  relatively 
coarse  mesh  employed  (Grid  G2),  the  LOC  boundary  method  incurs  significant  errors  (Fig.  4.20).  This 
situation  is  exacerbated  since  the  source  center  lies  on  one  of  the  boundaries  and  part  of  the  acoustic  energy 
is  dissipated  before  propagation.  The  problem  with  the  LOC  approach  is  overcome  through  optimization 
of  the  filter  coefficient  near  the  boundary.  As  Fig.  4.20  indicates,  the  higher-order  results  are  essentially 
recovered  by  the  optimized  LOC  (or  OLOC)  method  when  a/  is  set  to  0.49  for  the  second-  and  fourth-order 
filter  operators  (points  ‘2’  and  ‘3’).  These  results  reiterate  the  superior  performance  of  compact-based  filters 
over  standard  explicit  filters  since  the  latter  contain  no  mechanism  for  control  of  the  spectral  response.  This 
issue  is  of  considerable  importance  not  only  for  accuracy  as  already  shown,  but  also  for  robustness  (stability) 
of  the  overall  high-order  approach.  As  noted  in  reference  to  Fig.  3.2,  the  spectral  response  of  the  high-order 
boundary  filters  exhibit,  due  to  the  lack  of  symmetry,  undesirable  amplification  over  a  small  portion  of 
the  modified  wave  number  range.  For  a  fixed  order  of  accuracy,  this  behavior  is  more  pronounced  for  the 
standard  explicit  filters,  but  can  be  lessened  for  the  compact  filter  through  an  increase  in  ay.  A  practical 
example  where  this  control  is  beneficial  is  illustrated  in  Fig.  4.21  for  the  scattering  cylinder  configuration, 
in  which  results  with  the  C6F10  scheme  and  with  a  minimum  6th-order  boundary  filter  are  shown  near 
the  cylinder  surface.  The  computation  in  Fig.  4.21(a)  employed  a/  =  0.45  and  provided  a  stable  solution. 
However,  when  an  explicit  filter  is  used  (i.e.  a/  =  0),  spurious  oscillations  appear  near  the  surface  and 
ultimately  cause  numerical  instability.  By  increasing  ay  to  0.45  at  the  near-boundary  points  only,  these 
unwanted  oscillations  were  avoided.  This  suggests  that  the  observed  difficulties  are  due  to  the  asymmetric 
near-boundary  high-order  explicit  filter. 

4.3.4  Scattering  Multi-Domain  Application 

In  order  to  demonstrate  the  capability  of  the  present  numerical  approach  to  handle  multiple-domain  appli¬ 
cations,  the  previous  scattering  case  was  simulated  using  the  two-zone  overlap  configuration  shown  schemat¬ 
ically  in  Fig.  4.22.  The  original  single  grid  ( G2 ,  Table  4.2)  was  split  along  9  =  90°,  and  extra  £-lines  were 
added  to  form  a  five-point  overlap  (see  Fig.  4.1).  The  solutions  were  advanced  separately  on  each  domain, 
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Figure  4.20:  Effect  of  filter  boundary  treatment  on  directivity  at  r/D  =  5.  HO=Higher  order  boundary 
formulas,  LOC=Lower-order  centered,  OLOC=Optimized  lower-order  centered 


Figure  4.21:  Comparison  of  (a)  implicit  and  (b)  explicit  near-boundary  high-order  filters 
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Figure  4.22:  Schematic  of  mesh  overlap  for  multi-domain  cylinder  scattering  problem 


and  information  was  exchanged  at  the  overlap  points  in  the  manner  previously  discussed  in  Section  3.7. 
Solutions  were  obtained  using  the  sixth-order  compact  scheme  for  interior  points  along  with  fourth-  and 
fifth-order  compact  operators  at  the  boundary  and  next-to-boundary  points  respectively.  In  the  interior  of 
each  domain,  a  lOth-order  filter  was  employed.  Near  the  boundary,  two  different  choices  were  made,  the  first 
with  a  minimum  sixth-order  and  the  second  with  a  minimum  eighth-order  formula.  For  all  filter  operators, 
the  coefficient  a/  =  0.45  was  specified. 

Figure  4.23  displays  instantaneous  pressure  contours  in  the  vicinity  of  the  vertical  overlap.  It  is  apparent 
that  the  pressure  waves  cross  the  grid  interface  without  producing  any  noticeable  disruptions  of  the  interfer¬ 
ence  pattern.  It  should  be  noted  that  pressure  waves  generated  by  the  source  propagate  through  the  overlap 
region  in  an  oblique  direction  to  the  zonal  interface.  A  quantitative  comparison  between  the  single-  and 
multiple-domain  solutions  is  given  in  Fig.  4.24  in  terms  of  the  directivity  of  the  radiated  sound  at  r/D  =  5. 
The  directivity  obtained  with  either  6th-  or  8th-order  near-boundary  filter  formulas  in  the  overlap  region  is 
essentially  the  same  as  the  corresponding  single-domain  baseline  solution.  These  results  highlight  the  po¬ 
tential  of  the  present  high-order  methodology  for  domain-decomposition  applications  on  parallel  computers. 


4.3.5  Spherical  Pulse  on  a  3-D  Curvilinear  Mesh 


The  final  validation  case  considers  the  propagation  of  a  3-D  spherical  pulse  [53]  in  a  curvilinear  mesh.  An 
initial  pressure  pulse  is  prescribed  by 


lno(*2+*2+*2) 
V  =  Voo  +  e-e  a 


(4.3) 


where  e  =  0.01.  In  order  to  examine  metric  cancellation  errors,  a  three-dimensinal  curvilinear  mesh  is 
generated  using  the  following  equations: 


Vi,j,k  = 
%i,j,k  — 


min  +  Axa  [(i  -  1)+ 
Axsmnxy7rl'{-1'>Avo  sinrta;~7r(^~1)A~° 

J^y  J-*z 

ymin  +  A Vo  [(j  - 1)+ 
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IX  —  1  y  JL  —  1  KL  —  1 


(4.4) 
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Figure  4.23:  Pressure  contours  for  multi-domain  scattering  simulation  in  vicinity  of  mesh  overlap  region 
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Figure  4.24:  Comparison  of  directivity  at  r / D  =  5  for  multi-domain  and  single-domain  calculations 

with  the  specified  parameters  IL  =  JL  =  KL  =  61,  Ax  =  Ay  =  Az  =  1,  Lx  =  Ly  =  Lz  =  60  and 
nxy  =  nyz  =  ...  =  8.  These  parameters  yield  a  mesh  in  which  the  metric  identities  (e.g.  Eqn.  3.4)  are 
not  trivially  satisfied.  A  constant  £  plane  of  the  mesh  at  the  location  of  maximum  deformation  is  shown  in 
Fig.  4.25(a).  The  imposed  grid  undulations  are  resolved  with  approximately  15  PPW. 

The  pulse  propagation  problem  is  computed  with  a  At  =  0.004  using  the  CJ^FIO  algorithm  with  af  =  0.49. 
Figure  4.25(b)  displays  the  calculated  pressure  contours  on  a  plane  through  the  center  of  the  spherical  pulse 
at  t  =  10  for  the  case  when  the  metrics  are  evaluated  in  the  standard  fashion  of  Equation  3.3.  It  is 
apparent  that  significant  distortion  of  the  wave  front  occurs  due  to  the  lack  of  freestream  preservation  (i.e. 
metric  cancellation  errors).  The  perturbation  pressure  along  the  grid  line  i  =  j  =  31,  (Fig.  4.26)  indicates 
unacceptable  departure  from  the  theoretical  solution.  The  results  obtained  with  the  new  metric  evaluation 
procedure  of  Equation  3.5  exhibit  no  distortions  of  the  spherical  front  (Fig.  4.25(c))  and  are  in  excellent 
agreement  with  the  exact  solution  (Fig.  4.26). 

Calculations  with  C6F10  (not  shown)  displayed  reduced  sensitivity  to  the  choice  of  metric  evaluation 
procedure.  This  is  in  agreement  with  the  results  of  Ref.  [12],  where  metric  cancellation  errors  were  shown 
to  decrease  with  increasing  order  of  accuracy.  Nonetheless,  without  the  incorporation  of  Equation  3.5  in  the 
calculation  of  the  metric  relations,  all  solutions  on  this  highly  distorted  mesh  were  found  to  be  of  poor  quality. 
Even  for  relatively  benign  3-D  curvilinear  grids,  freestream  preservation  errors  may  swamp  acoustic  pressure 
levels  unless  proper  metric  evaluation  procedures  are  employed.  It  should  be  noted  that  analytic  metric 
evaluation  (even  if  available)  does  not  remedy  the  situation  [11,  1].  The  present  results  clearly  demonstrate 
that  the  superior  performance  of  the  high-order  method  can  be  extended  to  the  realm  of  general  curvilinear 
grids  making  the  approach  suitable  for  complex  configurations  including  moving/deforming  meshes  for  which 
good  grid  quality  cannot  usually  be  maintained. 
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Figure  4.25:  Three-dimensional  spherical  pulse,  (a)  Mesh,  (b)  and  (c)  Pressure  contours  with  standard  and 
new  metrics  respectively 


30 


Figure  4.26:  Effect  of  metric  evaluation  on  computed  pressure  along  line  through  spherical  pulse  at  t  =  10. 
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Chapter  5 


PREPROCESSING  TOOLS  FOR 
HIGH-ORDER  OVERSET-GRIDS 

5.1  Introduction 

A  new  family  of  preprocessing  tools  have  been  developed  for  use  with  the  high-orcler  overset-grid  (HO-OG) 
algorithm.  Some  of  these  tools  are  general  in  nature  and  have  been  implemented  to  improve  the  flexibility 
of  the  algorithm  as  currently  implemented,  while  others  are  specifically  needed  due  to  the  requirements  of 
the  high-order  method.  The  function  of  these  tools  is  to  quickly,  accurately  and  robustly  take  overset-grid 
systems  and  develop  the  necessary  input  for  use  with  the  high-order,  parallel  research  code  FDL3DI  [40]. 
The  approach  utilized  here  is  conceptually  similar  to  that  employed  by  the  BREAKUP  code  [54],  which 
some  important  modifications  to  address  unique  features  of  the  high-order  algorithm.  Section  5.2  addresses 
the  specific  preprocessing  elements  that  have  been  developed.  Section  5.3  details  the  flow  solver  and  specific 
modifications  that  have  been  made  to  it  to  accept  the  input  from  the  preprocessing  stage  as  well  as  to  improve 
its  overall  flexibility.  Finally,  two  examples  of  the  use  of  the  HO-OG  algorithm,  completed  with  the  newly 
developed  preprocessing  tools,  are  presented  in  Section  5.4.  These  include  a  two-dimensional  aeroacous- 
tic  application  involving  scattering  from  three  circular  cylinders  and  a  three-dimensional  electromagnetics 
problem  involving  scattering  from  a  single  sphere. 


5.2  Preprocessing  Elements 

The  approach  used  here  for  performing  HO-OG  simulations  is  shown  in  the  flowchart  in  Figure  5.1.  This 
section  will  focus  on  the  preprocessing  requirements,  while  the  following  section  will  address  issues  associated 
with  the  flow  solver  itself.  Postprocessing  will  not  be  specifically  addressed  in  this  work. 

For  the  purposes  of  the  discussion  here,  it  is  assumed  that  an  appropriate  overset-grid  system  [55]  has 
been  previously  generated  for  the  geometry  of  interest.  To  perform  a  computational  simulation  using  the 
overset-grid  method,  the  connectivity  between  overlapping  grids  must  first  be  established.  There  are  a  variety 
of  ways  to  do  this  [55],  but  here  this  task  is  performed  using  the  PEGASUS  Version  5  software  package.  [56] 
In  addition  to  the  grids  themselves  and  the  input  necessary  to  guide  the  operation  of  the  code,  PEGASUS 
requires  specification  of  the  boundary  conditions  that  will  be  enforced  in  the  solution  domain  by  the  flow 
solver.  Typically,  it  is  assumed  that  any  grid  boundaries  that  do  not  have  boundary  conditions  specified 
by  the  input  data  will  require  interpolation.  The  task  of  determining  the  appropriate  connectivity  between 
donor  and  receiver  grids  at  these  non-specified  boundaries  is  handled  by  PEGASUS. 

Upon  execution  of  PEGASUS,  a  file  is  generated  that  contains  the  connectivity  data  for  the  overset 
grid-level  topology.  This  data  is  contained  in  the  “XINTOUT”  file  and  includes  all  receiver  points  and  their 
corresponding  donor  points  and  offsets  (the  distance  between  the  base  point  of  the  donor  stencil  and  receiver 
point  in  the  computational  space  of  the  donor  grid) .  The  interpolation  offsets  are  calculated  using  an  inverse 
isoparametric  mapping  that  is  spatially  second-order  accurate  [57].  The  donor  point  is  actually  the  base  point 
identifying  the  computational  cell  containing  the  receiver  point.  The  eight  points  of  the  donor  cell  containing 
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Figure  5.1:  Workflow  diagram  for  the  high-order  overset-grid  process. 


the  receiver  point  provide  enough  data  to  perform  the  second-order  accurate  trilinear  interpolation  that  is 
typically  used  in  overset  methods.  [58]  A  second  required  file,  created  by  running  the  PEGASUS  utility 
pegplot,  will  be  referred  to  as  the  “XIBLANK”  file.  This  file  is  a  multi-block,  iblanked,  PLOT3D  grid  file 
where  the  iblank  value  for  each  receiver  point  is  equal  to  the  negative  value  of  the  corresponding  donor  grid. 

In  order  to  facilitate  the  application  of  the  overset-grid  approach  to  high-order  computational  simulations, 
several  additional  preprocessing  steps  are  necessary  once  PEGASUS  has  been  executed  and  the  XINTOUT 
and  XIBLANK  files  generated.  At  this  point  in  the  development  of  the  HO-OG  approach,  some  of  these  steps 
require  significant  user  intervention,  while  others  are  completely  automated.  The  main  preprocessing  code 
that  has  been  developed  to  perform  many  of  the  automated  tasks  in  the  HO-OG  process  is  referred  to  as 
BELLERO  (Bellerophon  was  the  mythological  Greek  hero  who  tamed  the  winged  horse  PEGASUS ,  upon 
whom  he  rode  to  slay  the  multi-component  monster  Chimera).  The  various  tasks  in  the  preprocessing  phase 
of  the  HO-OG  approach  will  now  be  discussed. 

5.2.1  Hole  Modification 

For  problems  requiring  PEGASUS  to  blank  out  points  in  one  or  more  grid  components,  the  first  step  is 
to  ensure  the  compatibility  of  the  resulting  holes  with  the  high-order  algorithm.  The  high-order  spatial 
algorithm,  and  in  particular  the  filtering  operation,  need  larger  stencils  than  are  typically  required  for  lower- 
order  simulations.  For  example,  a  2A"-order  interior  filter  of  the  type  employed  in  the  algorithm  require  a 
27V  +  1-point  stencil  [23].  Thus  a  10t,l-order  filter  requires  a  minimum  of  eleven  contiguous,  unblanked  points 
along  a  coordinate  direction  to  correctly  operate. 

The  requirement  to  maintain  adequate  stencil  sizes  for  the  high-order  differencing  and  filtering  operators 
places  some  restrictions  on  the  nature  of  the  holes  that  can  be  handled  with  the  HO-OG  approach.  The 
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Figure  5.2:  Situations  showing  linkage  between  holes  and  minimum  stencil  sizes 


obvious  one  is  that  any  holes  generated  must  not  result  in  stencils  with  fewer  than  some  some  minimum 
number  of  points  in  any  of  the  three  coordinate  directions.  This  minimum  value  is  dictated  by  the  order- 
of-accuracy  of  the  differencing  and/or  filtering  operator  to  be  used  by  the  solver.  Stencils  containing  points 
below  the  minimum  prescribed  value  may  arise  due  to  the  proximity  of  the  hole  to  a  exterior  domain  boundary 
(e.g.,  Figure  5.2(a)),  due  to  two  distinct  holes  in  close  proximity  to  each  other  (Figure  5.2(b)),  or  due  to  the 
presence  of  unblankecl  points  protruding  into  a  blanked  region  (Figure  5.2(c)).  This  later  situation  can  arise 
due  to  the  merging  of  holes  generated  by  two  or  more  distinct  bodies  in  close  proximity,  as  is  the  case  in 
Figure  5.2(c),  or  due  to  a  sufficiently  complex  grid  geometry  that  is  used  to  cut  the  hole. 

A  second  point  that  should  be  considered  when  working  with  holes  at  the  grid-topology  level  is  how 
the  holes  will  impact  the  partitioning  phase.  As  discussed  in  the  following  section,  the  minimum  stencil 
requirement  poses  a  constraint  on  how  the  partitioning  may  be  accomplished.  Thus,  the  hole  structure  at 
the  grid-topology  level  plays  an  important  role  in  the  quality  of  the  decomposition  at  the  block-topology 
level.  Here  it  will  suffice  to  make  the  general  statement  that  the  more  regular  the  hole  boundaries  are  (that 
is,  lacking  in  the  “stair-step”  behavior  where  the  hole  boundary  moves  in  multiple  coordinate  directions  at 
a  time),  the  more  flexibility  will  exist  when  decomposing  the  grids. 

The  situation  depicted  in  Figure  5.2(c)  arose  during  the  simulation  of  low-Reynolds  number  flow  over 
two  circular  cylinders  in  close  proximity.  [30]  Using  the  Level-2  interpolation  capability  [56]  of  PEGASUS , 
this  hole  pattern  was  cut  by  the  body-fitted  grids  surrounding  the  two  cylinders  into  the  background  grid. 
Identified  in  Figure  5.2(c)  are  some  of  the  regions  where  limited  stencil  support  existed.  The  HO-OG  method 
would  not  be  able  to  support  this  hole  structure,  and  thus  modifications  were  required  prior  to  performing 
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the  simulation.  The  modified  hole  structure  is  shown  in  Figure  5.2(d),  where  the  iblank  value  of  several 
points  has  been  modified  in  order  to  maintain  some  minimum  stencil  size  (fifteen  in  this  case). 

Currently,  the  modification  of  the  hole  points  must  be  done  in  manual  and  iterative  fashion.  First,  a 
utility  code  is  run  to  identify  and  report  the  locations  where  the  number  of  contiguous,  unblanked  points 
along  each  coordinate  direction  fall  below  some  user-specified  value.  A  second  utility  code  is  then  run  that 
allows  the  user  to  flip  the  values  of  the  iblank  array  in  the  XIBLANK  file  at  specified  points,  thus  potentially 
changing  the  shape  of  any  holes.  The  PEGASUS  utility  program  p3d2peg  is  then  used  to  enforce  the  new 
hole  map  in  PEGASUS,  which  is  run  again  to  establish  the  connectivity  for  this  modified  situation.  This 
manual  approach  of  check  stencils  /  modify  iblanks  /  rerun  PEGASUS  is  repeated  as  often  as  necessary 
in  order  to  obtain  a  grid  whose  connectivity  is  well-defined  (i.e.,  no  orphan  points  generated)  and  whose 
stencils  are  all  equal  to  or  larger  than  the  minimum  specified  size.  The  hole-cutting  parameters  in  PEGASUS 
may  also  be  modified  in  order  to  assist  in  generating  an  overset-grid  system  compatible  with  the  HO-OG 
minimum  stencil  requirement. 

Future  development  may  include  the  investigation  of  automated  methods  to  perform  the  necessary  hole 
modifications  with  minimal  user  input.  Such  an  approach  could  be  based  on  a  comparison  of  the  grid  points 
PEGASUS  would  prefer  to  blank  out  through  Level  2  interpolation  so  as  to  minimize  the  differences  in  cell 
size  and  orientation,  and  the  grid  points  that  PEGASUS  must  blank  out  due  to  their  location  inside  a  solid 
surface  boundary.  Such  a  comparison  would  identify  points  that  could  potentially  be  unblanked  in  order 
to  meet  HO-OG  minimum  stencil  requirements.  However,  this  approach  will  require  additional  research  to 
determine  its  feasibility  and  practicality  for  realistic  geometries. 

5.2.2  Automatic  Block  Partitioning 

The  parallelization  of  the  FDL3DI  code  was  performed  using  an  overset-grid  paradigm,  where  the  com¬ 
putational  grids  are  decomposed  into  overlapping  blocks  and  then  assigned  to  an  individual  processing 
element  [59].  While  the  decomposition  of  the  grids  into  blocks  is  relatively  straightforward  on  single-grid 
topologies,  there  are  load-balancing  issues  associated  with  domain  decomposition  approaches  on  general 
overset-grid  systems  [60] .  Also,  in  addition  to  being  satisfied  on  the  grid-level  topology,  the  minimum  stencil 
requirements  must  also  be  satisfied  at  the  block  level.  As  mentioned  in  the  previous  section,  the  presence 
of  holes  in  a  computational  domain  thus  introduces  additional  constraints  on  how  the  domain  may  be  par¬ 
titioned.  A  decomposition  that  is  both  load-balanced  and  meets  minimum  stencil  requirements  becomes 
difficult  to  obtain  manually  as  the  number  of  grids  and  the  complexity  of  the  holes  increases.  Thus,  au¬ 
tomated  methods  to  perform  the  decomposition  under  the  minimum  stencil  requirement  constraint  were 
investigated.  An  approach  was  developed  in-house  to  perform  this  task,  which  is  discussed  in  this  section. 

There  are  four  main  steps  to  the  automatic  partitioning  capability  developed  here.  The  first  step  is  to 
determine  the  number  of  blocks  into  which  each  grid  should  be  decomposed.  An  initial  estimate  of  the 
number  of  blocks  per  grid  is  performed  and  adjusted  until  the  total  number  of  blocks  matches  the  specified 
number  of  processors.  Next,  various  combinations  of  cuts  that  will  result  in  the  calculated  number  of  blocks 
on  each  grid  are  evaluated  to  determine  which  combination  produces  the  minimum  surface-to-volume  ratio 
in  order  to  minimize  communication  overhead.  This  is  done  through  a  simple  iterative  approach.  The  third 
step  is  to  then  apply  these  “optimal”  cuts  to  each  individual  grid  to  decompose  it. 

The  final  step  in  the  automated  decomposition  process  is  to  account  for  iblanking  and  enforce  minimum 
stencil  requirements.  This  is  done  by  first  flagging  every  point  in  the  grid  that  would  not  satisfy  the  minimum 
stencil  if  a  cut  were  to  be  made  through  it  in  a  given  direction  (thus  each  point  will  have  three  flag  values, 
one  for  each  of  the  three  coordinate  directions).  Thus,  the  shape  of  any  holes  will  have  an  impact  on  how  the 
grid  may  be  partitioned,  with  more  uniform  hole  boundaries  resulting  in  more  flexibility  when  decomposing 
the  grids.  Once  all  grid  points  have  been  flagged  and  the  individual  grids  have  been  cut,  a  check  is  performed 
to  find  all  cuts  which  cross  points  which  have  been  flagged  for  the  given  direction.  When  a  cut  is  found  to 
cross  a  point  which  violates  the  minimum  stencil  requirement,  the  cut  plane  is  moved  to  the  closest  valid 
plane.  After  all  of  the  necessary  cut  planes  have  been  moved,  the  remaining  planes  are  adjusted  in  order  to 
even  out  the  block  sizes  as  much  as  possible. 

An  alternate  way  to  account  for  iblanking  is  to  create  sub-grids.  When  creating  sub-grids  from  planes 
that  have  moved,  the  partitioning  process  is  identical  up  to  the  point  where  cut  planes  were  moved  in  order 
to  account  for  iblanking.  At  this  stage,  an  additional  step  is  added  to  create  sub-grids  which  satisfy  the 
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minimum  stencil.  These  sub-grids  are  then  partitioned  as  if  they  were  an  additional  grid  associated  with 
the  original  system  of  grids.  This  allows  for  highly  irregular  blocks  to  be  created,  which  otherwise  is  not 
possible  by  simply  making  cuts  which  extend  the  entire  plane  of  a  grid.  All  of  the  cut  planes  forced  to  move 
are  evaluated  as  potential  planes  to  create  sub-grids.  At  the  end  of  either  sub-grid  generation,  an  attempt 
is  made  to  even  out  the  size  of  the  sub-grids  where  possible.  This  is  done  by  shifting  cut  planes  where  the 
sub-grids  were  generated  on  the  sub-grid  which  contains  the  block  with  the  largest  volumes. 

The  above  process  is  still  in  the  developmental  phase.  It  has  been  successfully  demonstrated  for  a 
variety  of  straightforward  geometries  involving  holes  (see  Sec.  5.4),  but  needs  to  be  proven  for  more  complex 
configurations.  Additional  flexibility  in  allowing  the  user  to  specify  a  range  of  processors  to  consider  for 
obtaining  the  best  load-balanced  decomposition,  as  well  as  utilizing  other  established  domain  decomposition 
approaches  that  could  be  modified  to  account  for  iblanking,  are  concepts  under  consideration  for  future 
development . 

5.2.3  Block  Connectivity  and  Periodic  Boundary  Conditions 

Once  the  partitioning  for  the  grid  topology  has  been  specified,  the  block-level  connectivity  must  be  deter¬ 
mined.  There  are  three  components  to  this  step;  (1)  decomposing  the  grid-level  connectivity  data  provided 
by  PEGASUS  and  assigning  it  to  the  proper  blocks,  (2)  establishing  the  point-to-point  intragrid  connectivity 
requirements  such  that  the  blocks  decomposed  from  a  single  grid  can  communicate  with  each  other,  and  (3) 
establishing  the  point-to-point  connectivity  necessary  for  the  application  of  any  periodic  boundary  condi¬ 
tions.  Previously,  the  first  two  of  these  steps  were  handled  by  running  PEGASUS  on  the  decomposed  grid 
to  establish  the  necessary  block-level  connectivity.  However,  this  approach  had  several  drawbacks.  First, 
it  was  discovered  that  the  presence  of  block  boundaries  throughout  the  computational  domain  seemed  to 
hamper  the  ability  of  PEGASUS  to  create  high-quality  holes  using  its  automated  hole-cutting  process.  The 
work-around  developed  to  deal  with  this  issue  was  to  run  PEGASUS  to  perform  the  hole-cutting  on  the 
grid-level  topology,  and  then  rerun  it  using  the  resulting  hole-map  to  establish  the  block-level  connectiv¬ 
ity.  For  larger  grids  with  many  blocks,  the  time  required  to  rerun  PEGASUS  could  be  substantial  as  the 
connectivity  of  all  direct-injection  interpolation  points  had  to  be  established  using  the  same  algorithm  used 
for  the  donor/receiver  points  in  general  overlapping  regions  of  the  grids  (in  addition  to  the  reestablishing 
the  connectivity  of  the  already  calculated  overlaps  on  the  correct  block).  Also,  to  correctly  specify  the 
block-level  boundary  condition  input  data  required  by  PEGASUS  for  cases  with  increasingly  larger  num¬ 
ber  of  processors  required  significant  user  involvement  and  hence  was  time-consuming  and  prone  to  error. 
The  assignment  of  periodic  boundary  conditions  also  required  considerable  user  input,  with  the  necessary 
message-passing  constructs  hardwired  by  the  user  based  on  the  specifics  of  the  problem  and  the  particular 
decomposition  employed.  Finally,  the  resulting  connectivity  data  generated  through  this  process  was  valid 
only  for  a  particular  decomposition;  any  changes  (including  the  number  of  processors  to  run  the  case  on) 
required  a  significant  effort  on  the  part  of  the  user  in  recreating  the  PEGASUS  boundary  condition  data, 
rerunning  PEGASUS  on  the  new  block  topology,  and  modifying  the  code  to  correctly  handle  any  periodic 
boundary  conditions. 

Part  of  the  development  of  BELLERO  was  thus  geared  towards  finding  alternative  ways  to  determine 
block-level  connectivity.  This  issue  was  addressed  in  order  to  streamline  application  of  the  high-order  algo¬ 
rithm  to  new  configurations  of  interest  while  enhancing  the  flexibility  of  the  algorithm  by  relieving  the  user 
of  much  of  the  burden  associated  with  establishing  the  block-level  connectivity  for  a  particular  decomposi¬ 
tion.  One  of  the  key  features  of  BELLERO  is  its  ability  to  establish  the  block-level  connectivity  information 
for  a  specific  overset-grid  system  given  the  grid-level  connectivity  data  (obtained  from  PEGASUS)  and  an 
arbitrary  grid  partition  (either  user-supplied  or  determined  by  the  autopartitioning  capability  described  in 
the  previous  section).  The  point-to-point  connectivity  for  intragricl  boundaries  as  well  as  the  periodic  face 
boundaries  are  handled  in  a  similar  manner.  First,  the  linkages  between  periodic  faces  are  determined  using 
data  from  the  BCINP  namelist  in  the  PEGASUS  input  file.  Both  the  O-type  overlap  as  well  as  the  more 
general  case  where  the  periodic  faces  are  physically  separated  are  handled  by  BELLERO.  In  both  cases,  a 
boundary  condition  type  (variable  IBTYP  from  PEGASUS )  that  is  not  utilized  by  PEGASUS  is  selected 
to  flag  each  type  of  periodic  condition.  For  O-grid  periodicity,  each  pair  of  periodic  faces  is  flagged  using 
IBTYP=20  in  the  BCINP  namelist  input.  For  each  pair  of  non-O-type  periodic  faces,  a  unique  value  of  IBTYP 
>  200  is  selected  for  both  faces  to  provide  the  linkage  needed  to  establish  the  correct  donor/receiver  pairs. 
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DEFINE  GRIDS 


GRID-001  =  Body-fitted  cylindrical  grid  near  cylinder 
GRID-002  =  Cylindrical  near -background  grid 
GRID-003  =  Uake  grid 
GRID-004  =  Background  cartesian  grid 
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(a)  Example  BCINP  data  input  for  PEGASUS  with  periodic  (b)  Boundary  condition  output  generated  automatically 

boundary  conditions  specified  by  BELLERO  for  a  forty- four  block  decomposition  cor¬ 

responding  to  Figure  5.3(a) 

Figure  5.3:  Sample  boundary  condition  input  (to  PEGASUS )  and  output  (from  BELLERO ) 
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For  both  cases,  the  actual  regions  corresponding  to  the  periodic  overlaps  must  be  provided  by  the  user.  An 
example  of  BCINP  input  data  that  involves  the  use  of  both  types  of  periodic  boundary  conditions  is  given 
in  Figure  5.3(a).  This  data  corresponds  to  a  overset-grid  system  used  to  examine  the  transitional  flow  over 
a  circular  cylinder  at  Reo  =  3900  [31].  The  grid  system  consisted  of  four  grids;  a  body-fitted  cylindrical 
grid,  a  mid-field  intermediary  grid,  a  refined  wake  grid  and  a  background  cartesian  grid.  The  cylindrical  grid 
utilized  an  O-grid  periodic  face  on  the  Jmin  and  Jmax  faces  (IBDIR  ±1),  as  indicated  by  the  IBTYP  value 
of  20  on  each  of  these  faces.  In  the  L-direction  (spanwise),  periodic  boundary  conditions  were  employed, 
with  the  linkages  between  faces  indicated  by  the  various  IBTYP  values  greater  than  or  equal  to  200.  While 
in  this  case,  all  of  the  spanwise  periodic  linkages  connect  faces  on  the  same  grid,  the  IBTYP=2Cte  can  be 
used  to  link  together  faces  from  different  grids  as  required  by  the  topology  of  the  problem.  Conversely, 
the  IBTYP=20  (if  used)  must  be  specified  for  one  and  only  one  pair  of  faces  per  grid.  Also,  more  than  two 
faces  per  grid  can  be  specified  as  periodic  using  the  IBTYP=20x  flag  as  long  as  the  constraint  that  each  face 
matches  point-to-point  with  one  and  only  one  other  face  is  met.  In  the  case  where  two  or  more  adjacent 
faces  on  a  grid  are  periodic,  the  algorithm  will  enforce  all  periodic  conditions  simultaneously.  Thus,  in  the 
situation  on  a  computational  cube  where  each  face  is  periodic  with  its  opposite  face  (top-bottom,  left-right, 
front-back),  receiver  points  in  the  corners  of  the  cube  that  are  periodic  in  multiple  directions  will  be  directly 
paired  up  with  their  appropriate  donors  (top/right/front  with  bottom/left/rear,  etc.). 

Once  the  periodic  linkages  are  established,  the  algorithm  will  examine  all  of  the  block  faces  resulting 
from  the  partitioning  of  the  grid-level  topology.  If  a  block  face  is  either  an  intragrid  boundary  created  by 
the  partitioning  process,  or  if  its  an  exterior  grid  boundary  requiring  periodic  boundary  conditions,  then 
all  points  on  the  planes  that  are  to  receive  data  are  flagged  and  the  appropriate  donor  points  are  found. 
For  these  two  situations,  all  interpolation  is  direct-injection,  and  thus  the  offsets  are  set  identically  equal  to 
zero.  An  example  of  this  process  for  an  intragrid  boundary  is  shown  in  Figure  5.4,  which  shows  the  values 
in  the  iblank  array  on  an  original  grid  entity  (Figure  5.4(a))  and  then  the  iblank  values  on  one  of  the  block 
entities  (Block  20)  arising  from  the  decomposition  of  the  grid  (Figure  5.4(b)).  On  the  original  grid,  the 
iblank  values  in  the  region  shown  are  all  unity  indicating  that  these  points  are  field  points  as  determined  by 
PEGASUS.  After  finding  the  direct-injection  connectivity  between  blocks,  the  iblank  array  on  Block  20  has 
been  modified  to  show  which  blocks  now  provide  the  donor  information  along  the  block  boundaries,  indicated 
by  a  negative  iblank  value.  Thus,  on  the  portion  of  the  receiver  block  shown  in  Figure  5.4(b),  blocks  21,  24 
and  27  provide  the  necessary  donor  points.  As  shown  here,  a  five-point  overlap  region  with  a  two-point  fringe 
was  employed  between  blocks,  and  thus  a  set  of  non-communicating  overlap  points  [23]  exists  in  each  overlap 
region.  In  determining  whether  a  point  is  a  valid  donor,  the  intragrid  connectivity  algorithm  determines 
whether  it  is  also  a  receiver  point,  and  if  so,  eliminates  it  from  consideration  as  a  potential  donor.  Through 
this  process,  receiver  points  that  lie  in  regions  overlapped  by  more  that  two  blocks  (such  as  the  central  grey 
region  in  Figure  5.4(b))  will  receive  a  valid  donor  point  from  the  correct  neighboring  block.  For  points  where 
two  valid  direct-injection  donor  points  exist,  as  will  occur  along  edges  and  corners  of  blocks  at  NCO  points, 
the  algorithm  will  simply  select  the  point  on  the  lowest-numbered  block.  An  example  of  the  connectivity 
between  periodic  O-grid  faces  is  also  shown  in  Figs.  5.5,  which  demonstrates  the  connectivity  between  two 
blocks  (Blocks  4  and  8)  across  the  common  periodic  face. 

The  third  part  of  establishing  the  block-level  connectivity  involves  the  decomposition  of  the  grid-level  con¬ 
nectivity  provided  by  PEGASUS.  This  is  accomplished  by  finding  all  receiver  points  identified  by  PEGASUS 
that  reside  on  each  block,  and  then  determining  the  associated  donor  points  using  the  data  in  the  IBC  array 
from  PEGASUS.  Every  block  that  each  donor  point  could  reside  on  is  then  identified,  and  the  valid  donor 
point  that  results  in  the  largest  quality  factor  [61]  is  selected.  An  example  of  the  decomposition  of  grid-level 
connectivity  is  given  in  Figs.  5.6(a)  and  5.6(b).  The  original  values  in  the  iblank  array  associated  with  the 
receiver  grid  indicate  the  donor  grid  providing  the  data.  After  decomposition,  the  appropriate  block-level 
iblank  values  have  been  assigned.  The  donor/receiver  point  pairs  are  also  recast  into  the  block-level  indexing 
system  for  consistency  in  the  new  decomposed  topology. 

5.2.4  Boundary  Condition  Decomposition 

The  previous  approach  for  assigning  boundary  conditions  in  the  high-order  algorithm  also  required  consid¬ 
erable  user  input.  Boundary  conditions  were  assigned  on  the  block-level  topology  utilizing  logical  constructs 
hardwired  directly  into  the  code  to  assign  each  boundary  condition  to  the  appropriate  block.  While  satisfac- 
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(a)  Iblank  values  of  grid  prior  to  decomposition  (b)  Iblank  values  on  block  20  after  decomposition  showing 

the  donor  blocks  providing  data  at  the  block  boundaries 

Figure  5.4:  Example  of  intragrid  boundary  connectivity 


(a)  Iblank  values  of  block  4  along  O-grid  periodic  boundary 
(shared  with  block  8) 


(b)  Iblank  values  of  block  8  along  O-grid  periodic  boundary 
(shared  with  block  4) 


Figure  5.5:  Example  of  periodic  O-grid  boundary  connectivity 
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(a)  Iblank  values  of  grid  prior  to  decomposition 


(b)  Iblank  values  on  block  after  decomposition  showing 
the  decomposition  of  the  PEGASUS  intergrid  connectiv¬ 
ity  data 


Figure  5.6:  Example  of  intergrid  boundary  connectivity 


tory  for  basic  geometries  and/or  small  number  of  processors,  this  approach  was  clearly  self-limiting  as  more 
complex  configurations  were  considered  and  larger  numbers  of  processors  were  brought  to  bear  on  problems. 
Like  the  previous  approach  for  dealing  with  block-level  connectivity,  this  method  for  handling  boundary 
conditions  was  also  inflexible  in  the  sense  that  any  changes  in  the  partitioning  of  the  grids  would  require 
direct  modifications  to  the  code. 

To  increase  the  flexibility  of  the  algorithm  and  to  simplify  its  application  to  problems  involving  complex 
geometries  and/or  large  numbers  of  processors,  an  effort  was  made  relieve  the  user  of  the  burden  of  manually 
specifying  block-level  boundary  conditions.  As  PEGASUS  is  used  to  provide  grid-connectivity,  and  since 
PEGASUS  requires  boundary  condition  input  in  order  to  determine  which  grid  boundaries  will  receive 
interpolated  data,  it  was  decided  to  utilize  the  grid-level  boundary  condition  data  in  the  BCINP  namelist 
from  the  PEGASUS  input  file  as  the  starting  point  for  this  effort.  Consider  the  example  PEGASUS  BCINP 
namelist  input  as  given  in  Figure  5.3(a)  for  the  case  of  transitional  flow  over  a  circular  cylinder.  In  addition 
to  the  linkages  for  the  periodic  boundary  conditions,  all  other  grid-level  boundaries  that  require  boundary 
conditions  applied  in  the  flow  solver  are  specified  in  this  file.  The  new  boundary  condition  preprocessor  that 
is  part  of  BELLERO  now  takes  this  grid-level  boundary  condition  data  directly  from  the  PEGASUS  input 
file  and  determines  how  it  is  to  be  applied  at  the  block-level.  An  example  of  the  output  generated  by  the 
preprocessor  is  given  in  Figure  5.3(b)  for  a  forty-four  processor  decomposition.  This  data  is  then  used  by 
the  flow  solver  to  determine  how  to  apply  the  specified  boundary  conditions  to  the  particular  block-topology 
being  employed.  Note  that  the  boundary  conditions  associated  with  periodic  boundaries  in  Figure  5.3(b) 
are  not  represented  in  the  boundary  condition  decomposition  file  as  they  are  handled  through  the  domain 
connectivity  function. 

5.2.5  High-Order  Connectivity 

Once  it  has  determined  how  to  decompose  the  overset-grid  system  and  the  block-level  connectivity  and 
boundary  condition  assignments  have  been  generated,  BELLERO  then  begins  the  process  of  establishing 
high-order  connectivity  for  the  overset  grids  is  covered.  It  has  been  previously  shown  [62,  63]  that  the  use 
of  low-order  interpolation  is  not  sufficient  to  maintain  the  overall  accuracy  of  a  high-order  method.  Thus, 
extending  the  low-order  interpolation  data  provided  by  PEGASUS  to  higher  orders  is  an  important  part 
of  the  overall  HO-OG  algorithm.  Previous  work  [30]  determined  that  an  explicit  interpolation  procedure 
based  on  successive  one-dimensional  applications  of  Lagrangian  interpolation  formulae  was  well-suited  to 
maintain  overall  high-order  accuracy  within  the  HO-OG  framework,  and  that  is  the  only  interpolation 
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method  currently  employed  in  the  HO-OG  algorithm. 

The  establishment  of  high-order  connectivity  by  BELLERO  is  a  three-step  process.  In  the  first  step, 
the  second-order  stencils  generated  by  PEGASUS  are  expanding  such  that  the  dimension  of  the  stencil  in 
each  coordinate  direction  equals  the  specified  order  of  the  interpolation.  The  expansion  is  accomplished  such 
that  the  receiving  point  will  be  as  centered  as  possible  within  the  resulting  high-order  stencil.  The  second 
step  involves  the  calculation  of  the  high-order  interpolation  offsets  (the  mapping  of  the  physical  location  of 
the  receiver  point  in  the  computational  space  of  the  associated  donor  grid)  using  the  newly  created  high- 
order  stencils.  The  offsets  in  each  of  the  three  coordinate  directions  are  found  by  solving  three  isoparametric 
mapping  equations,  utilizing  the  second-order  offsets  found  by  PEGASUS  as  a  starting  point  for  the  iterative 
Newton’s  method  used  to  solve  the  equations.  The  final  step  is  the  calculation  of  the  Lagrangian  interpolation 
coefficients  associated  with  each  receiver  point  using  the  high-order  offsets.  A  total  of  3 N  coefficients  are 
required  to  fully  specify  a  three-dimensional,  Nth-OYder  interpolation  process  for  each  receiver  point.  More 
detail  on  the  determination  of  the  high-order  connectivity  can  be  found  in  Ref.  [30]. 

5.2.6  Hole  Management 

The  final  issue  addressed  by  BELLERO  is  that  of  hole  management.  The  concept  of  holes,  i.e.,  points  that 
are  excluded  from  consideration  by  the  flow  solver  due  to  their  location  inside  of  a  solid  surface  or  in  an 
overlapping  region  that  is  better  resolved  by  another  grid,  is  a  powerful  technique  in  overset-grid  methods. 
Its  utilization  greatly  enhances  the  overall  flexibility  of  the  approach  and  has  the  potential  to  simply  the 
generation  of  overset  grid  systems.  Hole  management  in  BELLERO  consists  of  three  elements;  the  identi¬ 
fication  of  hole  boundary  points  in  each  of  the  three  coordinate  directions,  the  insertion  of  the  appropriate 
one-sided  boundary  formulations  for  the  differencing  and  filtering  operations  at  all  points  impacted  by  the 
hole  boundary,  and  the  insertion  of  the  hole  decoupling  coefficients  at  all  points  within  the  hole  boundary. 
Each  of  these  three  elements  are  also  discussed  in  detail  in  Ref.  [30]. 

5.2.7  BELLERO  Input  and  Output 

To  summarize,  BELLERO  requires  as  input  four  files.  These  are  the  grid-level  XIBLANK  and  XINTOUT  files  that 
are  generated  by  PEGASUS  or  its  pegplot  utility,  a  PEGASUS- style  input  file  containing  the  BCINP  data  in 
order  to  process  periodic  boundary  conditions  and  decompose  other  boundary  conditions  for  the  block-level 
topology,  and  an  input  file  containing  input  data  to  manage  how  BELLERO  executes.  Upon  completion  of 
BELLERO ,  four  primary  files  are  generated  that  are  used  either  used  for  diagnostic  purposes  or  employed 
directly  into  the  flow  solver.  The  first  one  is  the  high-order,  block-level  XINTOUT  file.  This  file  is  structurally 
consistent  with  the  second-order,  grid-level  XINTOUT  file  generated  by  PEGASUS.  However,  it  is  based  on 
the  block-level  topology  and  thus  includes  all  of  the  necessary  point-to-point  communication  between  blocks 
as  well  as  any  periodic  boundary  conditions.  Also,  the  general  connectivity  data,  namely  the  base  donor 
points  and  interpolation  offsets,  have  been  upgraded  to  allow  for  high-order  interpolation.  The  second  file 
is  referred  to  as  the  “.X”  file  and  contains  some  additional  data  required  for  application  of  the  high-orcler 
interpolation  process,  such  as  the  size  of  the  stencil  and  the  Lagrangian  interpolation  coefficients  associated 
with  each  donor  point,  as  well  as  data  that  marks  the  start  and  finish  of  any  holes  along  each  constant 
coordinate  line  in  all  three  coordinate  directions.  The  third  file  is  the  boundary  condition  decomposition 
file,  an  example  of  which  was  shown  in  Figure  5.3(b).  This  file  allows  the  grid-level  boundary  conditions 
specified  in  the  PEGASUS  input  file  to  be  applied  on  the  block-level  topology  without  user  intervention. 
Finally,  the  fourth  file  is  the  decomposed  XIBLANK  file  containing  the  decomposed  grid  and  iblank  array  in 
PLOT3D  format.  This  file  is  very  useful  for  plotting  and  diagnostic  purposes. 


5.3  Flow  Solver 

The  baseline  flow  solver  for  use  with  BELLERO  is  FDL3DI  [40],  an  in-house  research  code  developed 
in  the  Computational  Sciences  Center  of  Excellence  of  the  Air  Force  Research  Laboratory  for  simulating 
the  Navier-Stokes  and  Euler  equations.  The  main  spatial  algorithm  is  based  on  fourth-  and  sixth-order 
compact  finite  differences  with  up  to  tenth-order  spatial  filters  to  maintain  stability  and  accuracy  by  remove 
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spurious  oscillations.  Available  time  integration  methods  include  both  an  explicit  Runge-Kutta  fourth- 
order  temporally  accurate  scheme  as  well  as  an  implicit,  approximately  factored  Beam- Warming  scheme 
of  up-to-third-orcler  temporal  accuracy  that  employs  diagonalization  and  Newton  subiterations  to  maintain 
efficiency  and  accuracy,  respectively.  The  code  has  been  equipped  with  a  general  overset-grid  capability 
using  Message  Passing  Interface  (MPI)  calls,  which  also  serves  as  a  domain  decomposition  paradigm  for 
utilization  on  parallel,  distributed  memory  platforms  [59] .  Additional  subroutines  have  been  incorporated  to 
perform  the  high-order  interpolation  using  the  coefficients  and  donor  stencils  calculated  by  BELLERO.  Also, 
hole  management  subroutines  have  been  included  to  read  in  BELLERO  generated  hole  data  and  insert  the 
appropriate  one-sided  differencing  and  filtering  formulas  at  hole  boundaries  as  well  as  decouple  the  points 
interior  of  the  hole.  Two  routines  were  added  to  process  the  new  boundary  condition  data  generated  by 
BELLERO ;  one  to  read  and  store  the  information  and  one  to  implement  the  boundary  conditions  in  a  general 
fashion.  Finally,  the  I/O  routines  were  modified  so  that  restart  data  could  be  read  in  using  the  block-level 
topology,  or  it  could  be  read  in  on  the  grid-level  topology  and  decomposed  internally  and  assigned  to  the 
appropriate  blocks. 

In  addition  to  this  code,  a  derivative  code  based  on  FDL3DI  was  developed  to  further  enhance  the 
flexibility  of  the  high-order  algorithm.  Titled  OHMS  (for  Overset  7/Igh-order  Maxwell  /Solver),  this  code 
employs  several  Fortran  90  features  such  as  allocatable  arrays  and  modular  code  structure  to  further  stream¬ 
line  problem  set-up.  All  subroutines  dependent  upon  the  equation  set  being  simulated  are  grouped  together 
in  modules,  which  enables  the  rapid  application  of  the  algorithm  to  a  variety  of  situations.  Three  equation 
sets  are  currently  using  this  approach;  the  Euler  equations  and  the  Maxwell  equations  (examples  of  which 
are  given  in  the  following  section),  and  the  Maxwell  equations  with  an  added  model  for  simulating  electro¬ 
magnetic  propagation  through  plasmas  [64].  Future  plans  include  incorporating  the  Navier-Stokes  equations 
into  the  OHMS  framework  as  well  as  adding  the  Beam- Warming  time-integration  scheme  (Runge-Kutta  is 
the  only  time  integration  scheme  implemented  in  the  OHMS  version  of  FDL3DI. 

5.4  Examples 

Two  validation  cases  are  examined  here  to  demonstrate  the  use  of  the  newly  developed  preprocessing  capa¬ 
bility  of  BELLERO  coupled  with  the  OHMS- based  version  of  FDL3DI.  Additional  validation  cases  may  be 
found  in  Ref.  [64]. 

5.4.1  Acoustic  Scattering  from  Three  Cylinders 

This  case  simulates  the  scattering  of  acoustic  waves  generated  by  a  spatially  distributed  source  from  three 
circular  cylinders.  It  was  one  of  the  benchmark  problems  in  the  Fourth  Computational  Aeroacoustic  Work¬ 
shop  [65],  and  serves  as  a  good  validation  case  due  to  the  presence  of  an  analytic  solution  [66]  for  comparison. 
This  particular  configuration  has  been  investigated  [32]  using  the  previous  approach  of  manually  decompos¬ 
ing  the  grids  and  boundary  conditions  down  to  the  block-level  as  well  as  rerunning  PEGASUS  after  the 
hole-cutting  phase  to  generate  all  block-level  connectivity.  Now,  PEGASUS  is  employed  to  cut  the  holes 
and  generate  the  grid-level  connectivity  only,  and  the  decomposition  of  the  grids  and  boundary  conditions 
and  the  generation  of  the  block-level  connectivity  is  handled  by  BELLERO. 

The  cylinder  configuration  shown  in  Figure  5.7(a)  is  gridded  using  three  meshes;  body-fitted  grids  about 
each  cylinder  embedded  in  a  background  Cartesian  mesh.  Additional  details  about  the  grid  system  employed 
as  well  as  the  parameters  of  the  acoustic  source  may  be  found  in  Ref.  [32].  Three  decompositions  were 
automatically  generated  using  BELLERO ,  including  sixteen,  thirty-two,  and  sixty-four  block  configurations. 
Pertinent  data  for  these  automatically  generated  partitions,  as  well  as  data  from  the  manually  generated, 
twenty-eight  block  partition  used  in  Ref.  [32],  is  shown  in  Table  5.1.  The  OHMS  code  with  the  Euler  equation 
module  was  employed,  using  sixth-order  interior  differencing,  tenth-order  spatial  filtering,  and  the  Runge- 
Kutta  time  integration  algorithm  with  At  =  0.002,  to  advance  the  solution  for  20,000  time  steps.  The  RMS 
fluctuating  pressure  profiles  on  the  surface  of  the  left  and  top  cylinders  were  collected  over  the  last  2,000 
time  steps  of  the  run  to  use  for  comparison  purposes.  The  results  from  the  BELLERO / OHMS  cases  are 
given  in  Figure  5.7(b)  and  show  a  good  agreement  between  the  various  cases,  the  previous  results,  and  the 
analytic  solution.  Note  that  generating  the  various  solutions  shown  here  only  required  that  the  parameter 
controlling  the  number  of  blocks  to  decompose  the  problem  into  be  modified  and  that  BELLERO  be  rerun. 
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(a)  Cylinder  configuration  (b)  RMS  fluctuating  pressure  profiles  on  surface 

of  cylinders 


(c)  Relative  speedup  for  various  block  decompositions 


Figure  5.7:  Results  from  acoustic  scattering  validation  problem 
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NP 

MBV 

LBC 

SVR 

MMR 

16 

65,510 

0.864 

0.698 

0.66 

32 

32,190 

0.875 

0.709 

0.67 

64 

14,976 

0.989 

0.727 

0.90 

28 

33,660 

0.967 

0.713 

0.64 

Table  5.1:  Block  decomposition  data  for  acoustic  scattering  case  (NP  =  number  of  processors,  MBV  = 
maximum  block  volume,  SVR  =  average  surface-to-volume  ratio,  MMR  =  minimum  volume  to  maximum 
volume  ratio) 


NP 

MBV 

LBC 

SVR 

MMR 

Wall-clock  time  (sec.) 

12 

195,112 

0.859 

0.107 

0.28 

9,700.3 

16 

171,820 

0.799 

0.082 

0.31 

13,381.8 

32 

86,469 

0.849 

0.118 

0.62 

8,015.3 

Table  5.2:  Block  decomposition  data  for  electromagnetic  scattering  case  (see  Table  5.1  for  column  headings) 


Conversely,  the  previous  methodology  would  have  required  generating  the  grid  indices  manually,  generating 
the  block-level  PEGASUS  boundary  condition  input,  rerunning  PEGASUS  to  regenerate  the  connectivity, 
modifying  the  code  parameters  to  correctly  dimension  the  arrays  and  the  block-level  logic  in  the  boundary 
condition  subroutine  and  recompiling  the  code,  and  generating  a  new  restart  file  for  the  block  configuration. 
This  process  would  have  to  be  repeated  for  each  block  topology  to  be  examined. 

While  no  formal  scalability  study  was  performed  at  this  time,  the  relative  speed-up  referenced  to  the 
sixteen  processor  case  on  ASC/MSRC  Compaq  ES45  platform  is  plotted  in  Figure  5.7(c).  When  the  par¬ 
titioning  is  done  automatically,  superlinear  speed-ups  are  obtained  through  this  range  of  processors.  Note 
the  poor  performance  of  the  manually  generated  partition  used  in  Ref.  [32]  compared  to  those  automatically 
by  BELLERO.  It  is  believed  that  this  is  due  to  the  large  variation  of  the  grid  dimensions  in  the  individual 
coordinate  directions  from  one  block  to  another  in  the  manual  partition,  which  could  be  causing  some  cache 
management  problems.  This  issue  is  currently  being  investigated. 

5.4.2  Electromagnetic  Scattering  from  Single  Sphere 

The  second  case  examined  here  is  the  scattering  of  electromagnetic  plane  waves  at  ka  =  4-7T  from  a  single 
sphere.  This  case  also  is  useful  as  a  validation  problem  due  to  the  existence  of  an  analytic  solution.  [67] 
The  sphere  was  gridded  using  the  overset  system  as  shown  in  Figure  5.4.2,  which  included  a  body-fitted 
polar  grid  with  the  poles  removed  around  the  sphere,  two  patch  grids  replacing  the  poles,  and  a  background 
Cartesian  grid.  BELLERO  was  used  to  handle  all  necessary  preprocessing  requirements  for  a  twelve,  sixteen 
and  thirty-two  block  configurations  whose  grid  parameters  are  given  in  Table  5.2.  The  FDL3DI- derived 
OHMS  code  was  used  as  the  flow  solver,  this  time  employing  the  Maxwell  equation  module.  The  same  solver 
configuration  was  employed  as  was  used  in  the  previous  aeroacoustic  case,  with  the  solution  being  advanced 
2000  time  steps  using  a  A t  =  0.005.  The  RMS  electric  field  on  the  surface  of  the  sphere  was  calculated  over 
the  last  1000  time  steps  for  each  decomposition,  and  the  results  are  shown  in  Figure  5.9  compared  to  the 
analytic  solution.  Again,  all  decompositions  match  the  analytic  solution  very  well.  The  wall-clock  times  for 
this  problem  are  also  included  in  Table  5.2.  For  this  case,  the  speed-ups  obtained  were  considerably  less 
than  the  previous  aeroacoustic  case,  and  in  fact  it  took  longer  to  run  this  case  on  sixteen  processors  than  it 
did  on  twelve  processors.  These  results  are  attributed  to  two  factors.  First,  the  irregular  hole  as  shown  in 
Figure  5.8(c)  limits  the  flexibility  of  the  decomposition  process  to  obtain  balanced  blocks  as  indicated  by  the 
parameters  in  Table  5.2.  Also,  the  sub-grid  approach  as  previously  discussed  for  automatically  decomposing 
the  grids  was  employed  for  the  sixteen  and  thirty-two  processor  cases.  This  resulted  in  blocks  whose  indices 
varied  considerably  from  block  to  block  in  the  three  coordinate  directions  (three  sample  block  sizes  from  the 
thirty-two  block  case  include  71  x  25  x  44,  19  x  83  x  44  and  41  x  19  x  83).  Thus,  cache  management  issues 
could  again  be  coming  into  play  in  these  cases. 
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(a)  Mid-plane  view 


(b)  Surface  grids 


Figure  5.8:  Overset  grid  system  used  for  scattering  from  sphere 


(a)  Comparison  of  scattered  RMS  electric  field  on  (b)  Scattered  RMS  electric  field  on  and  near 
surface  of  sphere  for  three  decompositions  sphere 


Figure  5.9:  Results  for  EM  scattering  from  sphere 
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Chapter  6 

COMMUNICATION  THROUGH 
PLASMA 


6.1  Introduction 

The  present  work  decouples  the  solution  of  the  non-equilibrium  re-entry  flow  from  the  problem  of  wave 
propagation.  A  second-order  Roe  solver[36]  is  utilized  to  obtain  the  overall  flow  conditions  and  then  the 
total  and  electron  number  densities  and  electron  temperature  are  used  to  solve  for  the  wave  propagation. 
Equations  for  conservation  of  electron  number  density  and  momentum  have  been  added  to  OHMS  (Overset 
High-Order  Maxwell  Solver)  developed  at  the  Air  Force  Research  Laboratory.  [64] 

6.2  Non-Dimensionalization  of  Governing  Equations 

The  equations  for  wave  propagation  through  an  ionized  medium  is  given  by  Maxwell’s  Equations  (in  the 
Rationalized  MKSA  system  [68]): 

e0c»tE  + J  -  V  x  (B/^o)  +  V  x  M  =  0  (6.1) 

dtB  +  V  x  E  =  0  (6.2) 

M  =  0  (6.3) 

where  the  current  is  given  by: 

J  =  <7snsvs  =  Jj  +  Je  (6-4) 

S 

with  Jj  =  <y(njVj  and  Je  =  — eneve.  The  number  density  equations  for  the  electrons,  ions  and  neutrals 
are[36,  69]: 

d t  (ne)  +  V  •  (neve)  =  ^  (wpej  +  u^)  (6.5) 

3 

dt  (rii)  +  V  •  (rijVj)  =  ij  +  d’ij)  (6-6) 

3 

dt  (n„)  +  V  •  (n„v„)  =  ^  (wpnj  +  w^)  (6.7) 

3 

where  A'A  is  the  production  of  species  s  from  species  j  and  is  the  destruction  of  species  s  into  species  j 
subject  to  the  constraint  that: 

H  +  ^sj)  =  0  (6-8) 

S  j 
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similarly,  the  species  momentum  equations  can  be  written  as [36,  69]: 


dt  (neve)  +  V  •  (neve  0  ve) 


nee 

(E  +  ve  x  B) 
me 


VPe 

me 


F  coll,e 

me 


+  +  V^ei) 

j 


(6.9) 


riiqi 

rtii 


^Pn  ^  Pcoll,i 
mn  mn 


E coll,i  (  -p,  •  d  \ 

m.  ■  2JV<<;  •  'V'-J 
j 

(6.10) 

-  +  (V^nj  +  VjW^j) 

(6.11) 

where  Fcou}S  is  a  force  term  for  species  s  due  to  collisions  and  the  last  term  on  the  right  hand  side  of 
Equations  (6.9)-(6.11)  is  the  transfer  of  momentum  due  to  creation  and  destruction  of  the  species. 

For  the  wave  propagation  problem,  the  variables  are  perturbed  about  the  steady  state  solution  obtained 
from  the  non-equilibrium  solution. 


E  =  E0  +  Ei 

B 

=  Bi 

—  ^e,0  +  ne,l 

Je 

Je,0 

Tli  =  0 

J, 

=  Ji,0 

=  ^n,0 

vn 

=  vn,0 

As  the  plasma  is  assumed  quasi-neutral,  n,,o  «  ne, o  =  no-  The  neutral  species  are  unaffected  by  the  wave 
propagation  due  to  the  lack  of  charge,  while  the  ions  can  be  considered  “frozen”  in  relation  to  the  electrons 
due  the  electron  mass  being  significantly  smaller  than  the  ion  mass.  For  the  same  reason,  the  collision  force 
term  is  only  a  function  of  electron  velocity: 


E  collie 

me 


-veneve 


where  ve  is  taken  as  an  average  collision  frequency  between  the  electrons  and  neutrals.  Electron-ion  collisions 
are  neglected  with  respect  to  the  electron-neutral  collisions  due  to  the  weak  ionization  of  the  flow  [70].  When 
the  steady-state  solution  is  subtracted  from  (6.1)-(6.11),  the  following  dimensional  equations  hold: 


eo<9*Ei  +  Je>i  —  V  x  (Bi/^o)  —  0 
9tBi  +  V  x  Ei  =  0 

dtneA  -  V  Je’1  =  0 


dt  Je>1  -  V 


e2ne 


E1 - (JeXBi) 

me 


7e-RTeVne,i  T 

6  e,l 


(6.12) 

(6.13) 

(6.14) 

(6.15) 


where  the  term  proportional  to  ne,iEo  is  neglected  due  to  the  quasi-neutral  approximation  of  the  plasma. 
Because  the  electron  plasma  waves  are  driven  by  the  Lorentz  force,  the  compressions  tend  to  be  anisotropic 
and  act  locally  one-dimensional,  so  qe  «  3  is  normally  assumed  [70,  71]. 

The  following  non-dimensionalizations  may  be  made: 


B*  =  Bi/6,  E*  =  Ei/£,  J*  =  J  e/J  (6.16) 

t*  =  t/r,  x*  =  x/L  (6-17) 

v:=ve/V  (6.18) 

n*  =  ne/n0 

,max  (6.19) 

to*  =  1  (6.20) 
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v*e  =  VeT 

where  the  non-dimensional  parameters  are  chosen  such  that: 


(6.21) 


L/t  =  c0  =  1/V^oeo 
B  =  £/cq 

~  £q£/t  C77.0, max'll 

which  reduces  to  the  non-dimensional  form: 


(6.22) 

(6.23) 

(6.24) 


dr  3 


* 

e,l 


at.Ej  +  j 


* 

e,l 


-  V*  x  BJ  =  0 


<9rB*  +  V*  x  E*  =  0 
drnh-jSV*-  Jh=0 


(co2pr2)  «EJ  -  /3(J:  x  B())  +  ^fTe*VX,i  - 


(6.25) 

(6.26) 

(6.27) 

(6.28) 


where  (3  =  V/co  is  a  characteristic  relativistic  Mach  number,  =  -\/ notmaxe2/ (eo me)  is  the  electron  plasma 
frequency  based  on  maximum  electron  density  and  vrms  =  \J?jKOe/rne  is  the  molecular  velocity  of  electrons 
with  temperature  scale  6e  (T*  =  Te/9e). 

Defining  new  nondimensional  variables:  f lp  =  ujpt  and  A 4th  =  V /vrrns ,  the  equation  for  the  electron 
current  may  be  rewritten  as: 


dt*j*eA-pv* 


=  nl  «Et-/3(j:  xBJ))  + 


M2th 


le,l 


*  T 

—  V(1A 


e,l 


(6.29) 


For  non-relativistic  electrons  one  may  assume  (3  «  1  and  the  non-linear  term  J*  (g>  J*  may  be  neglected. 
Since  the  current  simulation  assumes  that  there  is  no  externally  applied  B  field,  the  J*x  B|(  term  may  also 
be  neglected. 

Dropping  the  *’s  for  convenience,  the  equation  set  may  be  reduced  to: 


<9*Ei  +  J e,i  —  V  x  Bi  —  0 

(6.30) 

c^Bi  +  V  x  Ei  =  0 

(6.31) 

dtne,i  -/?V  Je,i  =  0 

(6.32) 

1  —  QpTleEi  1  +  2  ^^e,l  ^e^e,! 

(6.33) 

Note  that  the  divergence  of  current  in  Equation  (6.32)  cannot  be  dropped  until  M.2h  is  defined. 


6.3  Hypersonic  Model 

For  the  current  simulation  the  background  electrons  are  determined  using  the  same  second  order  Roe  solver 
developed  by  Josyula  and  Bailey.  [36]  The  hypersonic  code  utilizes  a  seven  species  model  consisting  of  N2, 
02,  N,  O,  NO,  NO+  and  e_.  The  mass  conservation  assumed  only  binary  collisions  and  the  mobility  was 
determined  using  a  constant  Lewis  number.  The  ions  experienced  ambipolar  diffusion  while  the  electrons 
were  assumed  to  be  in  Boltzmann  equilibrium  with  the  ions.  In  the  original  work,  the  wall  was  assumed 
non-catalytic  and  held  at  a  fixed  temperature  of  1500K;  however  in  the  present  work,  catalytic  walls  for  the 
electrons  was  assumed.  For  more  details  of  this  solution  methodology,  the  reader  is  directed  to  reference  [36]. 

The  free-stream  concentrations  in  the  original  code  were  too  high  for  the  atmospheric  conditions  at  61 
km.  While  the  values  were  still  small  enough  to  not  affect  the  hypersonic  shock  or  boundary  layer,  the 
electron  density  was  still  too  high  to  allow  waves  to  propagate  in  the  free-stream.  Figure  6.1  demonstrates 
this  insensitivity  to  the  free-stream  values  to  the  boundary  layer  along  the  stagnation  line. 
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x/Rn 


Figure  6.1:  Comparison  of  electron  density  along  the  stagnation  line  with  original  and  modified  free-stream 
values.  Distance  is  normalized  by  nose  radius  (Rn= 15.24  cm). 


Another  modification  to  the  code  was  the  addition  of  catalytic  boundary  conditions  at  the  wall.  The 
accumulated  space  charge  on  the  re-entry  vehicle  is  not  considered,  as  it  would  likely  require  not  only  solving 
the  plasma  sheath  in  more  detail,  but  also  the  trajectory  of  the  vehicle  to  account  for  its  history.  While 
the  exact  catalytic  nature  of  the  wall  is  not  known  in  general,  in  the  absence  of  space  charge  accumulation, 
the  surface  should  probably  be  fully  catalytic  to  at  least  the  electrons  and  ions [72].  Although  the  code 
solves  the  non-equilibrium  equations,  it  was  not  trivial  to  use  these  routines  to  determine  the  equilibrium 
conditions  needed  for  the  catalytic  surface.  Instead,  a  simple  algebraic  model  for  the  equilibrium  conditions 
was  incorporated. 

The  equilibrium  reactions  for  the  7  species  model  are  (cf.  [73]): 


with  equilibrium  constants  satisfying: 


N2^  N  +  N 

(6.34) 

O2  O  +  O 

(6.35) 

N  +  O^NO 

(6.36) 

N  +  O  ^  NO+  +  e~ 

(6.37) 

=  KN{T) 

Pn2 

(6.38) 

{P°)2  =  K0{T) 

Po2 

(6.39) 

=  Kno(T) 

PnPo 

(6.40) 

=  K  {T) 

PnPo 

(6.41) 

l  by: 

Cs  =  Ps/P 

(6.42) 

where  ps  is  the  partial  pressure  of  species  s,  and  P  is  the  total  pressure;  the  mole  fraction  for  atomic  nitrogen 
and  oxygen  are  related  as: 
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CN  =  \J  CN2  Kn  (T)/P  (6.43) 

C0  =  V Co2Kq(T)/P  (6.44) 

Since  at  1500  K,  there  is  very  little  dissociation  of  oxygen  and  nitrogen,  Cn2  ,  Co2  ~  const,  and  the  atomic 
species  can  be  approximated  as: 


Similarly,  for  the  mole  fraction  of  NO: 


Cn 


Co 


CN,Prcf 


Co,prBf 


CnO  =  CNO,PrBfCNCo 
and  utilizing  the  conservation  of  charge: 

Ce-  =  CNO+  =  CNO+,pref 


P 


P 


ref 


C’ref 

P 


(6.45) 

(6.46) 


(6.47) 


(6.48) 


To  ensure  that  the  total  number  of  species  is  conserved,  the  mole  fractions  for  molecular  nitrogen  and  oxygen 
are  calculated  from: 


Cn2  +  Cq2  +  Cno  +  Co  +  Cn  +  Cno+  +  Ce-  —  1 


(6.49) 


2Co2  +  Cno  +  Co  +  Cno+ 
2Cn2  +  Cno  +  Cn  +  Cno+ 


ftot 


(6.50) 


where  ftot  is  the  fraction  of  the  total  number  of  oxygen  atoms  divided  by  the  total  number  of  nitrogen  atoms. 

The  algebraic  approximations  used  the  equilibrium  values  for  T  =  1500° K  and  Pref  =  1  atm.  Equilibrium 
constants  were  calculated  from  values  in  the  JANAF  tables.  [74]  Figure  6.2  shows  the  agreement  between  the 
approximation  and  and  exact  solution  of  the  equilibrium  equations  along  the  body.  The  maximum  deviation 
from  the  exact  solution  is  less  than  0.12%. 


6.4  Results 

The  configuration  modeled  is  that  of  the  RAM  C-II  flight  at  an  altitude  of  61  km  and  Mach  number  of 
23.9.  [34,  35,  36,  64]  The  maximum  electron  temperature  in  the  flow  field  was  just  under  10000  K  and  the 
maximum  electron  density  was  approximately  1.4  x  1014/cm3.  The  length  scale  ( L )  was  taken  to  be  1  cm, 
the  temperature  scale  0e  to  be  11600K  and  the  electric  field  scale  was  chosen  to  be  equal  to  that  of  a  horn 
transmitting  at  1  kW  of  power.  Figure  6.3  shows  the  plasma  shielding  an  incoming  wave  transmitting  at 
9.21  GHz. 

The  non-equilibrium  solver  was  re-run  with  fully  catalytic  walls  and  walls  that  are  only  catalytic  to 
charged  particles.  Unlike  the  change  in  the  free-stream  conditions [64] ,  adding  a  catalytic  wall  boundary 
significantly  alters  the  electron  concentration  throughout  the  boundary  layer  as  seen  in  figure  6.4. 

The  horn  antenna  at  a  location  23.16  cm  along  the  axis  (location  2)  [34]  is  shown  in  figure  6.5.  While  the 
original  horn  was  circular,  an  equivalent  rectangular  horn  was  used  to  simplify  the  geometry.  The  vehicle 
grid  and  sub-domain  including  the  horn  can  be  seen  in  figure  6.6. 

The  presence  of  the  ionized  gas  blocking  signal  can  be  seen  in  figure  6.7.  The  main  factor  affecting  the 
signal  propagation  is  the  frequency  of  the  wave.  The  effect  of  the  catalytic  boundary  condition  was  also 
not  apparently  a  large  factor  as  seen  in  the  9.2  and  10.5  GHz  signals  in  figure  6.8.  The  difference  between 
that  and  an  equivalent  radiating  aperture,  while  quantitatively  different  closer  to  the  vehicle,  the  two  signals 
become  more  similar  further  away.  An  early  study[64]  indicated  that  the  solution  did  not  appear  to  be 
converged.  It  was  later  determined  that  in  the  three  grids,  the  antenna  apertures  were  slightly  different 
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Figure  6.2:  Comparison  of  atomic  oxygen  on  surface  of  vehicle  for  fully  catalytic  walls. 


Figure  6.3:  Plasma  shielding  incoming  9  GHz  wave.  Left:  no  plasma,  right:  plasma 
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Figure  6.4:  Comparison  of  electron  density  profile  normal  to  antenna  location  with  catalytic  and  non- 
catalytic  boundary  conditions. 


Figure  6.5:  RAM  C-II  re-entry  vehicle  with  locations  of  various  antennas  (ref:  NASA  TN  D-6062). 
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Amplitude 


Figure  6.6:  Computational  grids.  Left:  RAM  vehicle  prior  to  decomposition.  Right:  Grid  subset  near 
location  #2  for  embedded  horn. 


Figure  6.7:  Horn  driven  at  9.2  GHz  from  antenna  location  #2.  Left:  no  plasma,  right:  plasma 


Figure  6.8:  Different  frequencies  and  different  configurations  for  radiation  from  location  #2. 
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Figure  6.9:  Magnitude  of  instantaneous  electric  field  normal  to  the  center  of  the  radiating  antenna  for 
different  time-steps  and  grid  densities. 


sizes.  Figure  6.9  shows  the  grid  and  timestep  convergence  along  the  centerline  of  the  antenna  with  the 
non-catalytic  wall. 

It  turns  out  that  retaining  the  extra  terms  in  the  equations  had  a  negligible  effect  on  the  solution 
(difference  was  on  the  order  of  10~°  of  maximum  current).  In  hindsight,  this  should  not  have  been  surprising. 
For  the  current  configuration  with  the  horn  transmitting  at  1  kW  and  an  electron  temperature  ( KTe )  of 
leV,  f }p  ft:  22,  j^2~  ft:  8  and  ve  ss  0.01.  While  it  appears  that  the  electron  density  might  have  an  effect, 
the  factor  of  [3  in  (6.32)  is  small  enough  here  to  preclude  this.  Nor  would  increasing  the  horn  power  (which 
would  increase  (3)  help  in  retaining  the  extra  terms  in  this  linear  analysis.  This  is  because  the  relevant 
characteristic  parameter  is  which  is  independent  of  Mach  number  and  is  approximately  2  x  10~3. 

To  see  this,  an  alternate  perturbation  may  be  defined  for  the  number  density: 


ne  =  n6)  o  +  /?ne?i 

(6.51) 

in  which  case  (6.32)  and  (6.33)  simplify  to: 

1  —  V  •  Je,l  =  0 

(6.52) 

^t^e,  1  =  +  0thTeVfieil  ^eJe,l 

(6.53) 

where  [3th  =  Vrms/co-  Therefore,  in  order  to  see  an  effect  from  the  electron  pressure  term,  the  electron 
temperature  must  be  significantly  higher  than  is  generated  in  the  current  conditions. 
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Chapter  7 


CONCLUSIONS 


A  higher-order  compact-difference  based  method  has  been  adapted  to  simulate  aeroacoustic  phenomena  using 
the  Euler  and  Maxwell  equations.  The  components  of  the  scheme  include  fourth  and  sixth-order  Pade-type 
spatial  differencing  formulas  coupled  with  up  to  tenth-order  implicit  filters  which  are  required  to  maintain 
stability  in  the  presence  of  mesh  stretching  and  complex  boundary  conditions.  Significant  grid  stretching, 
common  in  practical  calculations,  is  shown  to  have  profound  impact  on  the  solution,  manifested  in  the 
generation  of  high-frequency  oscillations.  This  characteristic  is  employed  to  positive  effect  in  conjunction 
with  the  filtering  procedure  to  enhance  the  decay  of  acoustic  radiation  outside  the  domain  of  interest. 
The  potential  of  the  procedure  for  parallel  implementation  is  demonstrated  by  successful  application  to 
multidomain  scattering  cases.  Special  metric  procedures  are  shown  to  be  essential  in  the  simulation  of 
acoustic  phenomena  in  highly  curvilinear  three-dimensional  meshes. 

A  new  preprocessing  code  BELLERO  has  been  developed  to  automate  many  of  the  tasks  associated  with 
domain  decomposition  for  the  parallel,  high-order  overset-grid  flow  solver  FDL3DI.  Highlighted  capabilities 
of  BELLERO  include;  (1)  automatic  generation  of  the  grid  indices  for  the  domain  decomposition  taking 
into  account  minimum  stencil  requirements  for  the  high-order  algorithm,  (2)  automatic  generation  of  the 
block-level  connectivity  including  periodic  boundary  conditions,  (3)  automated  decomposition  of  the  grid- 
level  boundary  conditions  thus  eliminating  the  need  to  manually  specify  block-level  boundaries  in  the  code, 
and  (4)  calculation  of  high-order  interpolation  coefficients  and  management  of  hole  points  to  fully  implement 
the  high-orcler  overset-grid  (HO-OG)  approach.  Improvements  have  also  been  made  to  the  FDL3DI  solver 
itself  in  order  to  further  enhance  the  overall  flexibility  of  the  HO-OG  implementation.  The  new  capability 
was  demonstrated  for  benchmark  problems  of  acoustic  scattering  from  three  circular  cylinders  and  electro¬ 
magnetic  scattering  from  a  single  sphere.  Development  of  the  approach  is  continuing,  including  application 
on  more  complex  and  realistic  problems. 

The  ability  to  extend  the  wave  propagation  to  include  weakly  ionized  gases  has  been  demonstrated  in  a 
flexible  high-order  overset  framework.  It  has  been  shown  that  the  effect  of  adding  the  conservation  of  number 
density  to  the  equation  set  does  not  improve  the  solution  enough  to  justify  the  extra  time  and  storage  for 
re-entry  blackout  problems.  Surprisingly,  it  has  been  shown  that  changing  the  boundary  condition  to  make 
the  electrons  catalytic  at  the  wall  also  does  not  seem  to  have  a  very  large  effect. 
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NOMENCLATURE 


Uppercase  variables 


B 

B 

Cs 

E 

£ 

E 

F,G,H 

F  coll 

J 

J 

J 

Km 

L 

M 

Moo 

Mth 

N 

N2 

NO 

NO+ 

O 

02 

Ps 

R 

S 

T 

Te 

U 

u,v,w 

V 


Magnetic  Induction 
Reference  Magnetic  Field 
Mole-fraction  of  species  s 
Electric  Field  Intensity 
Reference  Electric  Field 
Energy 
Fluxes 

Force  due  to  collisions 

Current  Density 

Reference  current 

Jacobian  of  metric  transformation 

Equilibrium  constant  of  molecule  M 

Reference  length 

Magnetization 

Reference  Mach  number 

Thermal  Mach  number 

Atomic  Nitrogen 

Nitrogen  molecule 

Nitric  Oxide 

Nitric  Oxide  ion 

Atomic  Oxygen 

Oxygen  molecule 

Pressure  of  species  s 

Residual 

Source  vector 

Temperature 

Electron  temperature 

Solution  vector 

Velocity  components  in  generalized  coordinates 
Reference  velocity 


Lowercase  variables 

Co  Speed  of  light  in  free  space 

e  Electron,  or  charge  of  electron 

k  Wavenumber 

ns  Number  density  of  species  s 

ms  Mass  of  species  s 

p  Static  pressure 

qs  Charge  of  species  s 

t  Time 
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u,  v,  w  Velocity  components  in  Cartesean  coordinates 

vs  velocity  of  species  s 

Vrms  Electron  molecular  velocity 

ibgj  Destruction  of  species  s  into  species  j 

wU  Production  of  species  s  from  species  j 

x  A  general  coordinate  system 

x ,  y,  z  Cartesean  coordinates 


Greek  variables 

a  Tuning  parameter  for  compact-differences 

otf  Tuning  parameter  for  filters 

(3  Relativistic  Mach  number 

eo  Permittivity  of  free  space 

A  Change  of  variable  (i.e.  At  =  tn+1  —  tn) 

Po  Permeability  of  free  space 

vs  Collision  frequency  of  species  s 

p  Density 

t  Reference  time 

9e  Reference  temperature 

(f>  Generic  scalar 

77,  C  Generalized  curvilinear  coordinates 

ojp  Plasma  frequency 

Op  Non-dimensional  plasma  frequency 


Abbreviations 


AIAA 

ASC/MSRC 

BC 

C4 

C6 

CnFm 

CAA 

DRP 

E2 

ENO 

F8 

F 10 

FDTD 

HO 

HO-OG 

inc 

JANAF 

LOC 

MU  SC  L 

MVG 

OLOC 

NASA 

RAM 

RK 

RK4 

scat 

SFZ 


American  Institute  of  Aeronautics  and  Astronautics 

Aeronautical  Sciences  Center  Major  Shared  Resource  Center 

Boundary  Condition 

Fourth-order  compact-difference  scheme 

Sixth-order  compact-difference  scheme 

Compact-difference  of  order  n  with  filter  of  order  m 

Computational  Aeroacoustics 

Dispersion- Relation-Preserving 

Second-order  explicit  difference  scheme 

Essentially  Non-Oscillatory 

Eighth-order  implicit  filter 

Tenth-order  implicit  filter 

Finite-Difference  Time-Domain  (staggered  grid) 

High-Order 

High-Order  Overset-Grid 
Incident  field 

Joint  Army,  Navy,  Air  Force 
Lower-Order  Centered 

Monotone  Upstream-centered  Schemes  for  Conservation  Laws 

Method  of  minimization  of  group  velocity  errors 

Optimized  Lower-Order  Centered 

National  Aeronautics  and  Space  Administration 

Radio  Attenuation  Measurement 

Runge-Kutta  integration 

Standard  4^-order  Runge-Kutta 

Scattered  field 

Scattered  Field  Zone 
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TFZ 

TN 

tot 


Total  Field  Zone 
Technical  Note 
Total  field 


Codes  mentioned 


BELLERO 

BREAKUP 

FDL3DI 

OHMS 

PEGASUS 


Preprossing  tool  for  high-order  overset-grids  (developed  in-house) 

Code  to  decompose  overset  grids  developed  at  Sandia  National  Laboratories 
High-orcler  overset  code  to  solve  fluid  and  aeroacoustics  problems  (developed  in-house) 
Overset  High-order  Maxwell  Solver  derived  from  FDL3DI  (developed  in-house) 

Hole  cutting  program  for  overset  grids  developed  by  NASA 
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